most

MOST MATPOWER Optimal Scheduling Tool
   MDO = MOST(MDI)
   MDO = MOST(MDI, MPOPT)

   Solves a multiperiod, stochastic, contingency constrained, optimal
   power flow problem with linear constraints and unit commitment.
   Depending on inputs it may include DC power flow constraints or
   a simple total power balance condition.

   Inputs:
       MDI   MOST data structure, input
           (see MOST User's Manual for details)
       MPOPT   MATPOWER options struct, relevant fields are (default
               value in parens):
           verbose - see 'help mpoption'
           <solver specific options> - e.g. cplex, gurobi, etc,
                     see 'help mpoption'
           most.build_model (1) - build the MIQP, both constraints and
                   standard costs (not coordination cost) and store in
                   QP field of MDO
           most.solve_model (1) - solve the MIQP; if coordination
                   cost exists, update it; requires either 'most.build_model'
                   set to 1 or MDI.QP must contain previously built model
           most.resolve_new_cost (0) - use when MIQP is already built and
                   unchanged except for new coordination cost
           most.dc_model (1) - use DC flow network model as opposed to simple
                   generation = demand constraint
           most.fixed_res (-1) - include fixed zonal reserve contstraints,
                   -1 = if present, 1 = always include, 0 = never include
           most.q_coordination (0) - create Qg variables for reactive power
                   coordination
           most.security_constraints (-1) - include contingency contstraints,
                   -1 = if present, 1 = always include, 0 = never include
           most.storage.terminal_target (-1) - constrain the expected terminal
                   storage to target value, if present (1 = always, 0 = never)
           most.storage.cyclic (0) - if 1, then initial storage is a variable
                   constrained to = final expected storage; can't be
                   simultaneously true with most.storage.terminal_target
           most.uc.run (-1) - flag to indicate whether to perform unit
                   commitment; 0 = do NOT perform UC, 1 = DO perform UC,
                   -1 = perform UC if MDI.UC.CommitKey is present/non-empty
           most.uc.cyclic (0) - commitment restrictions (e.g. min up/down
                   times) roll over from end of horizon back to beginning
           most.alpha (0) - 0 = contingencies happen at beginning of period,
                   1 = at end of period
           most.solver ('DEFAULT') - see ALG argument to MIQPS_MATPOWER or
                   QPS_MATPOWER for details
           most.skip_prices (0) - skip price computation stage for mixed
                   integer problems, see 'help miqps_matpower' for details
           most.price_stage_warn_tol (1e-7) - tolerance on the objective fcn
               value and primal variable relative match required to avoid
               mis-match warning message, see 'help miqps_matpower' for details

   Outputs:
       MDO   MOST data structure, output
           (see MOST User's Manual for details)

0001 function mdo = most(mdi, mpopt)
0002 %MOST MATPOWER Optimal Scheduling Tool
0003 %   MDO = MOST(MDI)
0004 %   MDO = MOST(MDI, MPOPT)
0005 %
0006 %   Solves a multiperiod, stochastic, contingency constrained, optimal
0007 %   power flow problem with linear constraints and unit commitment.
0008 %   Depending on inputs it may include DC power flow constraints or
0009 %   a simple total power balance condition.
0010 %
0011 %   Inputs:
0012 %       MDI   MOST data structure, input
0013 %           (see MOST User's Manual for details)
0014 %       MPOPT   MATPOWER options struct, relevant fields are (default
0015 %               value in parens):
0016 %           verbose - see 'help mpoption'
0017 %           <solver specific options> - e.g. cplex, gurobi, etc,
0018 %                     see 'help mpoption'
0019 %           most.build_model (1) - build the MIQP, both constraints and
0020 %                   standard costs (not coordination cost) and store in
0021 %                   QP field of MDO
0022 %           most.solve_model (1) - solve the MIQP; if coordination
0023 %                   cost exists, update it; requires either 'most.build_model'
0024 %                   set to 1 or MDI.QP must contain previously built model
0025 %           most.resolve_new_cost (0) - use when MIQP is already built and
0026 %                   unchanged except for new coordination cost
0027 %           most.dc_model (1) - use DC flow network model as opposed to simple
0028 %                   generation = demand constraint
0029 %           most.fixed_res (-1) - include fixed zonal reserve contstraints,
0030 %                   -1 = if present, 1 = always include, 0 = never include
0031 %           most.q_coordination (0) - create Qg variables for reactive power
0032 %                   coordination
0033 %           most.security_constraints (-1) - include contingency contstraints,
0034 %                   -1 = if present, 1 = always include, 0 = never include
0035 %           most.storage.terminal_target (-1) - constrain the expected terminal
0036 %                   storage to target value, if present (1 = always, 0 = never)
0037 %           most.storage.cyclic (0) - if 1, then initial storage is a variable
0038 %                   constrained to = final expected storage; can't be
0039 %                   simultaneously true with most.storage.terminal_target
0040 %           most.uc.run (-1) - flag to indicate whether to perform unit
0041 %                   commitment; 0 = do NOT perform UC, 1 = DO perform UC,
0042 %                   -1 = perform UC if MDI.UC.CommitKey is present/non-empty
0043 %           most.uc.cyclic (0) - commitment restrictions (e.g. min up/down
0044 %                   times) roll over from end of horizon back to beginning
0045 %           most.alpha (0) - 0 = contingencies happen at beginning of period,
0046 %                   1 = at end of period
0047 %           most.solver ('DEFAULT') - see ALG argument to MIQPS_MATPOWER or
0048 %                   QPS_MATPOWER for details
0049 %           most.skip_prices (0) - skip price computation stage for mixed
0050 %                   integer problems, see 'help miqps_matpower' for details
0051 %           most.price_stage_warn_tol (1e-7) - tolerance on the objective fcn
0052 %               value and primal variable relative match required to avoid
0053 %               mis-match warning message, see 'help miqps_matpower' for details
0054 %
0055 %   Outputs:
0056 %       MDO   MOST data structure, output
0057 %           (see MOST User's Manual for details)
0058 
0059 
0060 %   MOST
0061 %   Copyright (c) 2010-2018, Power Systems Engineering Research Center (PSERC)
0062 %   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
0063 %   and Ray Zimmerman, PSERC Cornell
0064 %
0065 %   This file is part of MOST.
0066 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0067 %   See https://github.com/MATPOWER/most for more info.
0068 
0069 t0 = tic;
0070 
0071 %% default arguments
0072 if nargin < 2
0073     mpopt = mpoption;       %% use default options
0074 end
0075 
0076 verbose = mpopt.verbose;
0077 
0078 if verbose
0079     fprintf('\n=============================================================================\n');
0080     fprintf(  '          MATPOWER Optimal Scheduling Tool  --  MOST Version %s\n', mostver());
0081     fprintf(  '          A multiperiod stochastic secure OPF with unit commitment\n');
0082     fprintf(  '                       -----  Built on MATPOWER  -----\n');
0083     fprintf(  '  by Carlos E. Murillo-Sanchez, Universidad Nacional de Colombia--Manizales\n');
0084     fprintf(  '                  and Ray D. Zimmerman, Cornell University\n');
0085     fprintf(  '       (c) 2012-2019 Power Systems Engineering Research Center (PSERC)       \n');
0086     fprintf(  '=============================================================================\n');
0087 end
0088 
0089 %% if you want to do a normal solve, you have to create the QP
0090 if mpopt.most.solve_model && ~mpopt.most.resolve_new_cost
0091   mpopt = mpoption(mpopt, 'most.build_model', 1);
0092 end
0093 if ~mpopt.most.build_model && ~mpopt.most.solve_model
0094   error('most: Ah ... are you sure you want to do nothing? (either ''most.build_model'' or ''most.solve_model'' must be true)');
0095 end
0096 
0097 %% set up some variables we use throughout
0098 ng = size(mdi.mpc.gen, 1);
0099 nt = mdi.idx.nt;
0100 ns = length(mdi.Storage.UnitIdx);
0101 mdi.idx.ng = ng;
0102 mdi.idx.ns = ns;
0103 baseMVA = mdi.mpc.baseMVA;
0104 for t = 1:nt
0105   mdi.idx.nj(t) = length(mdi.tstep(t).OpCondSched);
0106 end
0107 if ~isfield(mdi, 'UC') || ~isfield(mdi.UC, 'CommitSched') || ...
0108         isempty(mdi.UC.CommitSched)
0109   if isfield(mdi, 'CommitSched') && ~isempty(mdi.CommitSched)
0110     warning('-----  most: MDI.CommitSched has moved to MDI.UC.CommitSched, please update your code.  -----');
0111     mdi.UC.CommitSched = mdi.CommitSched;
0112   else
0113     error('most: commitment schedule must be provided in MSPD_IN.UC.CommitSched');
0114   end
0115 end
0116 
0117 % set up model options
0118 UC = mpopt.most.uc.run;
0119 if UC == -1
0120   if isempty(mdi.UC.CommitKey)
0121     UC = 0;
0122   else
0123     UC = 1;
0124   end
0125 end
0126 mo = struct(...
0127     'DCMODEL',                      mpopt.most.dc_model, ...
0128     'IncludeFixedReserves',         mpopt.most.fixed_res, ...
0129     'SecurityConstrained',          mpopt.most.security_constraints, ...
0130     'QCoordination',                mpopt.most.q_coordination, ...
0131     'ForceCyclicStorage',           mpopt.most.storage.cyclic, ...
0132     'CyclicCommitment',             mpopt.most.uc.cyclic, ...
0133     'ForceExpectedTerminalStorage', mpopt.most.storage.terminal_target, ...
0134     'alpha',                        mpopt.most.alpha ...
0135 );
0136 if mo.IncludeFixedReserves == -1
0137   if isfield(mdi, 'FixedReserves') && isfield(mdi.FixedReserves(1,1,1), 'req')
0138     mo.IncludeFixedReserves = 1;
0139   else
0140     mo.IncludeFixedReserves = 0;
0141   end
0142 end
0143 if mo.SecurityConstrained == -1
0144   if isfield(mdi, 'cont') && isfield(mdi.cont(1,1), 'contab') && ...
0145           ~isempty(mdi.cont(1,1).contab)
0146     mo.SecurityConstrained = 1;
0147   else
0148     mo.SecurityConstrained = 0;
0149   end
0150 end
0151 if mo.ForceExpectedTerminalStorage == -1
0152   if isempty(mdi.Storage.ExpectedTerminalStorageAim) && ...
0153      isempty(mdi.Storage.ExpectedTerminalStorageMin) && ...
0154      isempty(mdi.Storage.ExpectedTerminalStorageMax)
0155     mo.ForceExpectedTerminalStorage = 0;
0156   else
0157     mo.ForceExpectedTerminalStorage = 1;
0158   end
0159 end
0160 if mo.IncludeFixedReserves && ~(isfield(mdi, 'FixedReserves') && ...
0161       isfield(mdi.FixedReserves(1,1,1), 'req'))
0162   error('most: MDI.FixedReserves(t,j,k) must be specified when MPOPT.most.fixed_res = 1');
0163 end
0164 if mo.SecurityConstrained && ~(isfield(mdi, 'cont') && ...
0165       isfield(mdi.cont(1,1), 'contab') && ~isempty(mdi.cont(1,1).contab))
0166   error('most: MDI.cont(t,j).contab cannot be empty when MPOPT.most.security_constraints = 1');
0167 end
0168 if mo.IncludeFixedReserves && mo.SecurityConstrained
0169   warning('most: Using MPOPT.most.fixed_res = 1 and MPOPT.most.security_constraints = 1 together is not recommended.');
0170 end
0171 if mo.ForceExpectedTerminalStorage == 1;
0172   if mo.ForceCyclicStorage
0173     error('most: storage model cannot be both cyclic and include a terminal target value; must change MPOPT.most.storage.cyclic or MPOPT.most.storage.terminal_target');
0174   end
0175   if ns && isempty(mdi.Storage.ExpectedTerminalStorageAim) && ...
0176            isempty(mdi.Storage.ExpectedTerminalStorageMin) && ...
0177            isempty(mdi.Storage.ExpectedTerminalStorageMax)
0178     error('most: MDI.Storage.ExpectedTerminalStorageAim|Min|Max cannot all be empty when MPOPT.most.storage.terminal_target = 1');
0179   end
0180   if ~isempty(mdi.Storage.ExpectedTerminalStorageAim)
0181     mdi.Storage.ExpectedTerminalStorageMin = mdi.Storage.ExpectedTerminalStorageAim;
0182     mdi.Storage.ExpectedTerminalStorageMax = mdi.Storage.ExpectedTerminalStorageAim;
0183   end
0184 end
0185 if UC && (~isfield(mdi.UC, 'CommitKey') || isempty(mdi.UC.CommitKey))
0186   error('most: cannot run unit commitment without specifying MDI.UC.CommitKey');
0187 end
0188 
0189 if ns
0190   if isempty(mdi.Storage.InitialStorage)
0191     error('most: Storage.InitialStorage must be specified');
0192   end
0193   if isempty(mdi.Storage.TerminalChargingPrice0)
0194     mdi.Storage.TerminalChargingPrice0 = mdi.Storage.TerminalStoragePrice;
0195   end
0196   if isempty(mdi.Storage.TerminalDischargingPrice0)
0197     mdi.Storage.TerminalDischargingPrice0 = mdi.Storage.TerminalStoragePrice;
0198   end
0199   if isempty(mdi.Storage.TerminalChargingPriceK)
0200     mdi.Storage.TerminalChargingPriceK = mdi.Storage.TerminalStoragePrice;
0201   end
0202   if isempty(mdi.Storage.TerminalDischargingPriceK)
0203     mdi.Storage.TerminalDischargingPriceK = mdi.Storage.TerminalStoragePrice;
0204   end
0205   if isempty(mdi.Storage.MinStorageLevel)
0206     error('most: Storage.MinStorageLevel must be specified');
0207   else
0208     MinStorageLevel = mdi.Storage.MinStorageLevel;
0209   end
0210   if size(MinStorageLevel, 1) == 1 && ns > 1    %% expand rows
0211     MinStorageLevel = ones(ns, 1) * MinStorageLevel;
0212   end
0213   if size(MinStorageLevel, 2) == 1 && nt > 1    %% expand cols
0214     MinStorageLevel = MinStorageLevel * ones(1, nt);
0215   end
0216   if isempty(mdi.Storage.MaxStorageLevel)
0217     error('most: Storage.MaxStorageLevel must be specified');
0218   else
0219     MaxStorageLevel = mdi.Storage.MaxStorageLevel;
0220   end
0221   if size(MaxStorageLevel, 1) == 1 && ns > 1    %% expand rows
0222     MaxStorageLevel = ones(ns, 1) * MaxStorageLevel;
0223   end
0224   if size(MaxStorageLevel, 2) == 1 && nt > 1    %% expand cols
0225     MaxStorageLevel = MaxStorageLevel * ones(1, nt);
0226   end
0227   if isempty(mdi.Storage.InEff)
0228     InEff = 1;                      %% no efficiency loss by default
0229   else
0230     InEff = mdi.Storage.InEff;
0231   end
0232   if size(InEff, 1) == 1 && ns > 1  %% expand rows
0233     InEff = ones(ns, 1) * InEff;
0234   end
0235   if size(InEff, 2) == 1 && nt > 1  %% expand cols
0236     InEff = InEff * ones(1, nt);
0237   end
0238   if isempty(mdi.Storage.OutEff)
0239     OutEff = 1;                     %% no efficiency loss by default
0240   else
0241     OutEff = mdi.Storage.OutEff;
0242   end
0243   if size(OutEff, 1) == 1 && ns > 1 %% expand rows
0244     OutEff = ones(ns, 1) * OutEff;
0245   end
0246   if size(OutEff, 2) == 1 && nt > 1 %% expand cols
0247     OutEff = OutEff * ones(1, nt);
0248   end
0249   if isempty(mdi.Storage.LossFactor)
0250     LossFactor = 0;                     %% no losses by default
0251   else
0252     LossFactor = mdi.Storage.LossFactor;
0253   end
0254   if size(LossFactor, 1) == 1 && ns > 1 %% expand rows
0255     LossFactor = ones(ns, 1) * LossFactor;
0256   end
0257   if size(LossFactor, 2) == 1 && nt > 1 %% expand cols
0258     LossFactor = LossFactor * ones(1, nt);
0259   end
0260   if isempty(mdi.Storage.rho)
0261     rho = 1;                        %% use worst case by default (for backward compatibility)
0262   else
0263     rho = mdi.Storage.rho;
0264   end
0265   if size(rho, 1) == 1 && ns > 1    %% expand rows
0266     rho = ones(ns, 1) * rho;
0267   end
0268   if size(rho, 2) == 1 && nt > 1    %% expand cols
0269     rho = rho * ones(1, nt);
0270   end
0271   if isempty(mdi.Storage.InitialStorageLowerBound)
0272     if mo.ForceCyclicStorage        %% lower bound for var s0, take from t=1
0273       mdi.Storage.InitialStorageLowerBound = MinStorageLevel(:, 1);
0274     else                            %% Sm(0), default = fixed param s0
0275       mdi.Storage.InitialStorageLowerBound = mdi.Storage.InitialStorage;
0276     end
0277   end
0278   if isempty(mdi.Storage.InitialStorageUpperBound)
0279     if mo.ForceCyclicStorage        %% upper bound for var s0, take from t=1
0280       mdi.Storage.InitialStorageUpperBound = MaxStorageLevel(:, 1);
0281     else                            %% Sp(0), default = fixed param s0
0282       mdi.Storage.InitialStorageUpperBound = mdi.Storage.InitialStorage;
0283     end
0284   end
0285   
0286   LossCoeff = mdi.Delta_T * LossFactor/2;
0287   beta1 = (1-LossCoeff) ./ (1+LossCoeff);
0288   beta2 = 1 ./ (1+LossCoeff);
0289   beta3 = 1 ./ (1/(1-mo.alpha) + LossCoeff);
0290   beta4 = mo.alpha/(1-mo.alpha) * beta2 .* beta3;
0291   beta5 = beta1 ./ beta2 .* (beta3 + beta4);
0292   beta2EtaIn      = mdi.Delta_T * beta2 .* InEff;
0293   beta2overEtaOut = mdi.Delta_T * beta2 ./ OutEff;
0294   beta3EtaIn      = mdi.Delta_T * beta3 .* InEff;
0295   beta3overEtaOut = mdi.Delta_T * beta3 ./ OutEff;
0296   beta4EtaIn      = mdi.Delta_T * beta4 .* InEff;
0297   beta4overEtaOut = mdi.Delta_T * beta4 ./ OutEff;
0298   diagBeta2EtaIn1      = spdiags(beta2EtaIn(:,1),      0, ns, ns);
0299   diagBeta2overEtaOut1 = spdiags(beta2overEtaOut(:,1), 0, ns, ns);
0300   diagBeta3EtaIn1      = spdiags(beta3EtaIn(:,1),      0, ns, ns);
0301   diagBeta3overEtaOut1 = spdiags(beta3overEtaOut(:,1), 0, ns, ns);
0302   diagBeta4EtaIn1      = spdiags(beta4EtaIn(:,1),      0, ns, ns);
0303   diagBeta4overEtaOut1 = spdiags(beta4overEtaOut(:,1), 0, ns, ns);
0304 end
0305 if ~isfield(mdi.idx, 'ntds') || isempty(mdi.idx.ntds) || ~mdi.idx.ntds
0306   ntds = 0;
0307   nzds = 0;
0308   nyds = 0;
0309 else
0310   ntds = mdi.idx.ntds;
0311   nzds = size(mdi.dstep(1).A, 1);
0312   nyds = size(mdi.dstep(1).D, 1);   % # of outputs of dynamical system
0313 end
0314 mdi.idx.ntds = ntds;
0315 mdi.idx.nzds = nzds;
0316 mdi.idx.nyds = nyds;
0317 
0318 %% define named indices into data matrices
0319 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0320     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0321 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0322     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0323     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0324 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0325     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0326     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0327 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0328 [CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
0329     CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
0330     CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
0331     CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
0332     CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
0333     CT_MODCOST_X] = idx_ct;
0334 
0335 %% check that bus numbers are equal to indices to bus (one set of bus numbers)
0336 nb  = size(mdi.mpc.bus, 1);
0337 if any(mdi.mpc.bus(:, BUS_I) ~= (1:nb)')
0338     error('most: buses must be numbered consecutively in bus matrix; use ext2int() to convert to internal ordering')
0339 end
0340 
0341 % Make data tables with full # of cols and add also pseudo OPF results to
0342 % be able to run printpf on them
0343 mdi.mpc.bus(:, MU_VMIN) = 0;
0344 mdi.mpc.gen(:, MU_QMIN) = 0;
0345 mdi.mpc.branch(:,MU_ANGMAX) = 0;
0346 mdi.mpc.f = 0;
0347 mdi.mpc.et = 0;
0348 mdi.mpc.success = 1;
0349 
0350 if mpopt.most.build_model
0351   if verbose
0352     fprintf('- Building indexing structures.\n');
0353   end
0354 
0355   %% save model options in data structure
0356   mdi.DCMODEL                               = mo.DCMODEL;
0357   mdi.IncludeFixedReserves                  = mo.IncludeFixedReserves;
0358   mdi.SecurityConstrained                   = mo.SecurityConstrained;
0359   mdi.QCoordination                         = mo.QCoordination;
0360   mdi.Storage.ForceCyclicStorage            = mo.ForceCyclicStorage;
0361   mdi.Storage.ForceExpectedTerminalStorage  = mo.ForceExpectedTerminalStorage;
0362   mdi.UC.run                                = UC;
0363   mdi.UC.CyclicCommitment                   = mo.CyclicCommitment;
0364   mdi.alpha                                 = mo.alpha;
0365   if ~isfield(mdi, 'OpenEnded'), mdi.OpenEnded = 1; end
0366 
0367   if UC
0368     % Make sure MinUp and MinDown are all >= 1
0369     if any(mdi.UC.MinUp < 1) && any(mdi.UC.MinUp < 1)
0370         error('most: UC.MinUp and UC.MinDown must all be >= 1');
0371     end
0372     % Unless something is forced off in mdi.CommitKey, or as a result of
0373     % not fulfilling its mdi.UC.MinDown in early periods, it should be available
0374     % for commitment and thus a contingency including its outage should not
0375     % be deleted.
0376     mdi.UC.CommitSched = (mdi.UC.CommitKey >= 0);   % Treat anything but -1 as on.
0377     if ~mdi.UC.CyclicCommitment
0378       for i = 1:ng
0379         if mdi.UC.InitialState(i) < 0
0380           nn = mdi.UC.MinDown(i) + mdi.UC.InitialState(i);  % time to go before startup
0381           if nn > 0
0382             mdi.UC.CommitSched(i, 1:nn) = 0;
0383           end
0384         elseif mdi.UC.InitialState(i) > 0
0385           nn = mdi.UC.MinUp(i) - mdi.UC.InitialState(i);    % time to go before shutdown
0386           if nn > 0
0387             mdi.UC.CommitSched(i, 1:nn) = 1;
0388           end
0389         end
0390       end
0391     end
0392   end
0393   % From now on, mdi.UC.CommitSched has zeros for stuff that is definitely
0394   % off, and ones for stuff that might be on, so those zeros can be used to
0395   % trim off contingencies that won't happen.
0396   % Start by creating the base flow data for all scenarios
0397   for t = 1:nt
0398     for j = 1:mdi.idx.nj(t)
0399       mpc = mdi.mpc;
0400       mpc.gen(:, GEN_STATUS) = mdi.UC.CommitSched(:, t);
0401 
0402       %% for backward compatibility, putting time dependent energy offer
0403       %% data in offer(t).gencost is deprecated, please use profiles
0404       if isfield(mdi.offer(t), 'gencost') && ~isempty(mdi.offer(t).gencost)
0405         mpc.gencost = mdi.offer(t).gencost;
0406       end
0407 
0408       if ~isempty(mdi.tstep(t).OpCondSched(j).tab)
0409         changelist = unique(mdi.tstep(t).OpCondSched(j).tab(:, CT_LABEL));
0410         for label = changelist'
0411           mpc = apply_changes(label, mpc, mdi.tstep(t).OpCondSched(j).tab);
0412         end
0413       end
0414       mdi.flow(t,j,1).mpc = mpc;
0415       mdi.idx.nb(t,j,1) = size(mdi.flow(t,j,1).mpc.bus, 1);
0416       mdi.idx.ny(t,j,1) = length(find(mdi.flow(t,j,1).mpc.gencost(:, MODEL) == PW_LINEAR));
0417     end
0418   end
0419   % Then continue to create contingent flow scenarios, deleting any
0420   % redundant contingencies (i.e., decommitting a gen or branch when its
0421   % status is guaranteed to be off). No rows are deleted from gen or branch,
0422   % but the number of contingencies can indeed change.
0423   for t = 1:nt
0424     % Set default ramp reserve mask, if not provided
0425     if ~isfield(mdi.tstep(t), 'TransMask') || isempty(mdi.tstep(t).TransMask)
0426       mdi.tstep(t).TransMask = ones(size(mdi.tstep(t).TransMat));
0427     end
0428     % First get current step's scenario probabilities
0429     if t == 1
0430       scenario_probs = mdi.tstep(1).TransMat; % the probability of the initial state is 1
0431     else
0432       scenario_probs = mdi.tstep(t).TransMat * mdi.CostWeights(1, 1:mdi.idx.nj(t-1), t-1)'; % otherwise compute from previous step base cases
0433     end
0434     mdi.StepProb(t) = sum(scenario_probs); % probability of making it to the t-th step
0435     if mdi.SecurityConstrained
0436       for j = 1:mdi.idx.nj(t)
0437         [tmp, ii] = sort(mdi.cont(t,j).contab(:, CT_LABEL)); %sort in ascending contingency label
0438         contab = mdi.cont(t,j).contab(ii, :);
0439         rowdecomlist = ones(size(contab,1), 1);
0440         for l = 1:size(contab, 1)
0441           if contab(l, CT_TABLE) == CT_TGEN  && contab(l, CT_COL) == GEN_STATUS ...
0442               && contab(l, CT_CHGTYPE) == CT_REP && contab(l, CT_NEWVAL) == 0 ... % gen turned off
0443               && mdi.flow(t,j,1).mpc.gen(contab(l, CT_ROW), GEN_STATUS) <= 0    % but it was off on input
0444            rowdecomlist(l) = 0;
0445           elseif contab(l, CT_TABLE) == CT_TBRCH && contab(l, CT_COL) == BR_STATUS ...
0446               && contab(l, CT_CHGTYPE) == CT_REP && contab(l, CT_NEWVAL) == 0 ... % branch taken out
0447               && mdi.flow(t,j,1).mpc.branch(contab(l, CT_ROW), BR_STATUS) <= 0  % but it was off on input
0448             rowdecomlist(l) = 0;
0449           end
0450         end
0451         contab = contab(rowdecomlist ~= 0, :);
0452         mdi.cont(t, j).contab = contab;
0453         clist = unique(contab(:, CT_LABEL));
0454         mdi.idx.nc(t, j) = length(clist);
0455         k = 2;
0456         for label = clist'
0457           mdi.flow(t, j, k).mpc = apply_changes(label, mdi.flow(t, j, 1).mpc, contab);
0458           ii = find( label == contab(:, CT_LABEL) );
0459           mdi.CostWeights(k, j, t) = contab(ii(1), CT_PROB);
0460           mdi.idx.nb(t, j, k) = size(mdi.flow(t, j, k).mpc.bus, 1);
0461           mdi.idx.ny(t, j, k) = length(find(mdi.flow(t, j, 1).mpc.gencost(:, MODEL) == PW_LINEAR));
0462           k = k + 1;
0463         end
0464         mdi.CostWeights(1, j, t) = 1 - sum(mdi.CostWeights(2:mdi.idx.nc(t,j)+1, j, t));
0465         mdi.CostWeights(1:mdi.idx.nc(t,j)+1, j, t) = scenario_probs(j) * mdi.CostWeights(1:mdi.idx.nc(t,j)+1, j, t);
0466       end
0467     else
0468       for j = 1:mdi.idx.nj(t)
0469         mdi.idx.nc(t, j) = 0;
0470         mdi.CostWeights(1, j, t) = scenario_probs(j);
0471       end
0472     end
0473   end
0474 
0475   % Compute adjusted (for alpha) cost weights for objective function
0476   if mdi.SecurityConstrained && mdi.alpha ~= 0
0477     for t = 1:nt
0478       for j = 1:mdi.idx.nj(t)
0479         mdi.CostWeightsAdj(1, j, t) = mdi.CostWeights(1, j, t);
0480         for k = 2:mdi.idx.nc(t,j)+1
0481           mdi.CostWeightsAdj(k, j, t) = (1-mdi.alpha) * mdi.CostWeights(k, j, t);
0482           mdi.CostWeightsAdj(1, j, t) = mdi.CostWeightsAdj(1, j, t) + mdi.alpha * mdi.CostWeights(k, j, t);
0483         end
0484       end
0485     end
0486   else
0487     mdi.CostWeightsAdj = mdi.CostWeights;
0488   end
0489 
0490   % If UC, also need to (possibly) modify gencosts so that each fm(p) at
0491   % p=0 is zero, so that fm(p) + u*c00 is equal to the original f(p). This
0492   % portion of the cost is valid at t if the unit is commited there, but
0493   % only for base scenarios and contingencies in which this particular unit
0494   % is not ousted, so must be careful later when probability-weighting the
0495   % corresponding u(i,t)!
0496   if UC
0497     if ~isfield(mdi.UC, 'c00') || isempty(mdi.UC.c00)   % if not empty assume
0498       mdi.UC.c00 = zeros(ng, nt);                       % contains correct info!
0499       for t = 1:nt
0500         for j = 1:mdi.idx.nj(t)
0501           for k = 1:mdi.idx.nc(t,j)+1
0502             mpc = mdi.flow(t,j,k).mpc;
0503             c00tjk = totcost(mpc.gencost, zeros(ng,1));
0504             mdi.UC.c00(:, t) = mdi.UC.c00(:, t) + mdi.CostWeightsAdj(k, j, t) * c00tjk;
0505             c0col = COST + mpc.gencost(:,NCOST) - 1;
0506             ipoly = find(mpc.gencost(:, MODEL) == POLYNOMIAL);
0507             ipwl  = find(mpc.gencost(:, MODEL) == PW_LINEAR);
0508             ii = sub2ind(size(mpc.gencost), ipoly, c0col(ipoly));
0509             mpc.gencost(ii) = mpc.gencost(ii) - c00tjk(ipoly);
0510             for i = ipwl'
0511               jj = COST+1:2:COST+2*mpc.gencost(i,NCOST)-1;
0512               mpc.gencost(i, jj) = mpc.gencost(i, jj) - c00tjk(i);
0513             end
0514             mpc.fixed_gencost = c00tjk;
0515             mdi.flow(t,j,k).mpc = mpc;
0516           end
0517         end
0518       end
0519     end
0520   end
0521   
0522   % Build variable indexing mechanism
0523   % Find total number of flows, buses and ny variables; (including offline gens)
0524   mdi.idx.nf_total = 0;
0525   mdi.idx.nb_total = 0;
0526   for t = 1:nt
0527     for j = 1:mdi.idx.nj(t)
0528       mdi.idx.nf_total = mdi.idx.nf_total + (1 + mdi.idx.nc(t,j));
0529       for k = 1:mdi.idx.nc(t,j)+1
0530         mdi.idx.nb_total = mdi.idx.nb_total + size(mdi.flow(t, j, k).mpc.bus, 1);
0531         ii = find(mdi.flow(t,j,k).mpc.gencost(:, MODEL) == PW_LINEAR);
0532         mdi.idx.ny(t,j,k) = length(ii);
0533       end
0534     end
0535   end
0536   mdi.idx.ns_total = ns * mdi.idx.nf_total;
0537   % Variable order resembles that of several C3SOPFs stacked together,
0538   % including the internally generated y variables, and then all of the
0539   % other new variables that are specific to HP, but excluding qg on one hand, and pc,
0540   % rp and rm since these now are common across several scenarios. So create first a matrix
0541   % of indices to the beginning of each c3sopf cell's vars.  Include the
0542   % mechanism for adding theta variables if we want to create DC flow restrictions.
0543   % Then start assigning the start and end indices for variables in each
0544   % c3sopf cell
0545   om = opt_model;
0546   nj_max = max(mdi.idx.nj);
0547   nc_max = max(max(mdi.idx.nc));
0548   Ing = speye(ng);
0549   Ins = speye(ns);
0550   if mdi.DCMODEL
0551     om.init_indexed_name('var', 'Va', {nt, nj_max, nc_max+1});
0552   end
0553   om.init_indexed_name('var', 'Pg', {nt, nj_max, nc_max+1});
0554   om.init_indexed_name('var', 'dPp', {nt, nj_max, nc_max+1});
0555   om.init_indexed_name('var', 'dPm', {nt, nj_max, nc_max+1});
0556   om.init_indexed_name('var', 'y', {nt, nj_max, nc_max+1});
0557   if mdi.IncludeFixedReserves
0558     om.init_indexed_name('var', 'R', {nt, nj_max, nc_max+1});
0559   end
0560   for t = 1:nt
0561     for j = 1:mdi.idx.nj(t)
0562       % first all angles if using DCMODEL
0563       for k = 1:mdi.idx.nc(t,j)+1
0564         if mdi.DCMODEL
0565           iref = find(mdi.flow(t,j,k).mpc.bus(:,BUS_TYPE) == REF);
0566           if verbose && length(iref) > 1
0567             errstr = ['\nmost: Warning: Multiple reference buses.\n', ...
0568               '           For a system with islands, a reference bus in each island\n', ...
0569               '           may help convergence, but in a fully connected system such\n', ...
0570               '           a situation is probably not reasonable.\n\n' ];
0571             fprintf(errstr);
0572           end
0573           Va0 = mdi.flow(t,j,k).mpc.bus(:,VA)*pi/180;
0574           Va_max = Inf(mdi.idx.nb(t,j,k), 1);
0575           Va_min = -Va_max;
0576           Va_min(iref) = mdi.flow(t,j,k).mpc.bus(iref,VA)*pi/180;
0577           Va_max(iref) = Va_min(iref);
0578           
0579           om.add_var('Va', {t,j,k}, mdi.idx.nb(t,j,k), Va0, Va_min, Va_max);
0580         end
0581       end
0582       % All active injections in c3sopf cell
0583       for k = 1:mdi.idx.nc(t,j)+1
0584         mpc = mdi.flow(t,j,k).mpc;
0585         genmask = mpc.gen(:,GEN_STATUS) > 0;
0586         p0 = genmask .* mpc.gen(:,PG) / baseMVA;
0587         if UC       % relax bounds here, enforced by uPmax, uPmin constraints
0588           pmin = genmask .* (min(mpc.gen(:, PMIN) / baseMVA, 0) - 1);
0589           pmax = genmask .* (max(mpc.gen(:, PMAX) / baseMVA, 0) + 1);
0590         else        % enforce bounds here, subject to flow's GEN_STATUS
0591           pmin = genmask .* mpc.gen(:, PMIN) / baseMVA;
0592           pmax = genmask .* mpc.gen(:, PMAX) / baseMVA;
0593         end
0594         om.add_var('Pg', {t,j,k}, ng, p0, pmin, pmax);
0595       end
0596       if mdi.IncludeFixedReserves
0597         for k = 1:mdi.idx.nc(t,j)+1
0598           mpc = mdi.flow(t,j,k).mpc;
0599           r = mdi.FixedReserves(t,j,k);
0600           nrz = size(r.req, 1); %% number of reserve zones
0601           if nrz > 1
0602               r.rgens = any(r.zones);   %% mask of gens available to provide reserves
0603           else
0604               r.rgens = r.zones;
0605           end
0606           r.igr = find(r.rgens);        %% indices of gens available to provide reserves
0607           ngr = length(r.igr);          %% number of gens available to provide reserves
0608           %% check data for consistent dimensions
0609           if size(r.zones, 1) ~= nrz
0610               error('most: the number of rows in FixedReserves(%d,%d,%d).req (%d) and FixedReserves(%d,%d,%d).zones (%d) must match', t, j, k, nrz, t, j, k, size(r.zones, 1));
0611           end
0612           if size(r.cost, 1) ~= ng && size(r.cost, 1) ~= ngr
0613               error('most: the number of rows in FixedReserves(%d,%d,%d).cost (%d) must equal the total number of generators (%d) or the number of generators able to provide reserves (%d)', t, j, k, size(r.cost, 1), ng, ngr);
0614           end
0615           if isfield(r, 'qty') && size(r.qty, 1) ~= size(r.cost, 1)
0616               error('most: FixedReserves(%d,%d,%d).cost (%d x 1) and FixedReserves(%d,%d,%d).qty (%d x 1) must be the same dimension', t, j, k, size(r.cost, 1), t, j, k, size(r.qty, 1));
0617           end
0618           %% convert both cost and qty from ngr x 1 to full ng x 1 vectors if necessary
0619           if size(r.cost, 1) < ng
0620               r.original.cost = r.cost;     %% save original
0621               cost = zeros(ng, 1);
0622               cost(r.igr) = r.cost;
0623               r.cost = cost;
0624               if isfield(r, 'qty')
0625                   r.original.qty = r.qty;   %% save original
0626                   qty = zeros(ng, 1);
0627                   qty(r.igr) = r.qty;
0628                   r.qty = qty;
0629               end
0630           end
0631           mdi.FixedReserves(t,j,k).rgens = r.rgens;
0632           mdi.FixedReserves(t,j,k).igr   = r.igr;
0633           if isfield(r, 'original')
0634               mdi.FixedReserves(t,j,k).original = r.original;
0635           end
0636           mdi.FixedReserves(t,j,k)       = r;   %% for cost & qty (now that fields match)
0637           Rmax = Inf(ngr, 1);               %% bound above by ...
0638           kk = find(mpc.gen(r.igr, RAMP_10));
0639           Rmax(kk) = mpc.gen(r.igr(kk), RAMP_10);   %% ... ramp rate and ...
0640           kk = find(r.qty(r.igr) < Rmax);
0641           Rmax(kk) = r.qty(r.igr(kk));      %% ... stated max reserve qty
0642           Rmax = Rmax / baseMVA;
0643           om.add_var('R', {t,j,k}, ngr, [], zeros(ngr, 1), Rmax);
0644         end
0645       end
0646       % All deltaP plus in c3sopf cell
0647       for k = 1:mdi.idx.nc(t,j)+1
0648         om.add_var('dPp', {t,j,k}, ng, [], zeros(ng,1), []);
0649       end
0650       % All deltaP minus in c3sopf cell
0651       for k = 1:mdi.idx.nc(t,j)+1
0652         om.add_var('dPm', {t,j,k}, ng, [], zeros(ng,1), []);
0653       end
0654       % All y variables in c3sopf cell - even if not committed.  There must
0655       % be a fixed cost associated with u(t,i,j) such that if u(t,i,j) = 0,
0656       % then the cost interpolated from the (x,y) pairs is zero, and if
0657       % u(t,i,j) = 1, then the fixed cost plus that interpolated from the
0658       % (x,y) pairs is as desired.
0659       for k = 1:mdi.idx.nc(t,j)+1
0660         om.add_var('y', {t,j,k}, mdi.idx.ny(t,j,k), [], [], []);
0661       end %
0662     end % for j
0663   end % for t
0664   % Continue with pc, rpp, rpm, one set for each time period
0665   om.init_indexed_name('var', 'Pc', {nt});
0666   om.init_indexed_name('var', 'Rpp', {nt});
0667   om.init_indexed_name('var', 'Rpm', {nt});
0668   for t = 1:nt
0669     om.add_var('Pc', {t}, ng);
0670     %% non-negativity on Rpp and Rpm is redundant, leave unbounded below
0671     %% (except where gen is off-line for all j and k)
0672     Rpmin = -Inf(ng,1);
0673     off = ones(ng,1);
0674     for j = 1:mdi.idx.nj(t);
0675       for k = 1:mdi.idx.nc(t,j)+1
0676         off = off & mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) <= 0;
0677       end
0678     end
0679     Rpmin(off == 1) = 0;
0680     om.add_var('Rpp', {t}, ng, [], Rpmin, mdi.offer(t).PositiveActiveReserveQuantity/baseMVA);
0681     om.add_var('Rpm', {t}, ng, [], Rpmin, mdi.offer(t).NegativeActiveReserveQuantity/baseMVA);
0682   end
0683   % Now load following ramping reserves.  In open ended problem, we need to
0684   % specify nt-1 ramping reserves, those needed to transition 1-2, 2-3, ..
0685   % (nt-1)-nt . The initial ramp constraint (from t=0 to t=1) is data, not a
0686   % variable.  But in terminal state (at t=nt+1) or cyclical problems (when
0687   % the t=nt to t=1 transition is also considered) we need nt ramping
0688   % reserves.
0689   if ~mdi.OpenEnded
0690     mdi.idx.ntramp = nt;
0691   else
0692     mdi.idx.ntramp = nt - 1;
0693   end
0694   om.init_indexed_name('var', 'Rrp', {mdi.idx.ntramp});
0695   om.init_indexed_name('var', 'Rrm', {mdi.idx.ntramp});
0696   for t = 1:mdi.idx.ntramp
0697     ramp30 = mdi.flow(t,1,1).mpc.gen(:,RAMP_30)*2*mdi.Delta_T;
0698     om.add_var('Rrp', {t}, ng, [], zeros(ng,1), ...
0699         min(mdi.offer(t).PositiveLoadFollowReserveQuantity, ramp30)/baseMVA);
0700   end
0701   for t = 1:mdi.idx.ntramp
0702     ramp30 = mdi.flow(t,1,1).mpc.gen(:,RAMP_30)*2*mdi.Delta_T;
0703     om.add_var('Rrm', {t}, ng, [], zeros(ng,1), ...
0704         min(mdi.offer(t).NegativeLoadFollowReserveQuantity, ramp30)/baseMVA);
0705   end
0706   % Continue with storage charge/discharge injections, one of each
0707   % for each flow; first all charge injections, then all discharge
0708   % injections
0709   om.init_indexed_name('var', 'Psc', {nt, nj_max, nc_max+1});
0710   om.init_indexed_name('var', 'Psd', {nt, nj_max, nc_max+1});
0711   for t = 1:nt
0712     for j = 1:mdi.idx.nj(t)
0713       for k = 1:mdi.idx.nc(t,j)+1
0714         if ns
0715           om.add_var('Psc', {t,j,k}, ns, [], [], zeros(ns,1));
0716         end
0717       end
0718     end
0719   end
0720   if ns
0721     for t = 1:nt
0722       for j = 1:mdi.idx.nj(t)
0723         for k = 1:mdi.idx.nc(t,j)+1
0724           om.add_var('Psd', {t,j,k}, ns, [], zeros(ns,1), []);
0725         end
0726       end
0727     end
0728   end
0729   % Continue with storage upper and lower bounds, one for each time period
0730   % and unit
0731   om.init_indexed_name('var', 'Sp', {nt});
0732   om.init_indexed_name('var', 'Sm', {nt});
0733   if ns
0734     for t = 1:nt
0735       om.add_var('Sp', {t}, ns, [], [], MaxStorageLevel(:,t)/baseMVA);
0736     end
0737     for t = 1:nt
0738       om.add_var('Sm', {t}, ns, [], MinStorageLevel(:,t)/baseMVA, []);
0739     end
0740   end
0741   % Possible initial storage quantities when using cyclic storage dispatch
0742   % so that initial storage = expected terminal storage is a constraint
0743   if ns && mdi.Storage.ForceCyclicStorage
0744     om.add_var('S0', ns, [], ...
0745         mdi.Storage.InitialStorageLowerBound / baseMVA, ...
0746         mdi.Storage.InitialStorageUpperBound / baseMVA);
0747   end
0748   % If there is a dynamical system with non-null state vector,
0749   % add those states here
0750   if nzds
0751     om.init_indexed_name('var', 'Z', {ntds});
0752     for t = 1:ntds
0753       if t == 1
0754         zmin = mdi.z1;
0755         zmax = mdi.z1;
0756       else
0757         zmin = mdi.dstep(t).zmin;
0758         zmax = mdi.dstep(t).zmax;
0759       end
0760       z0 = (zmax - zmin) / 2;
0761       om.add_var('Z', {t}, nzds, z0, zmin, zmax);
0762     end
0763   end
0764   % Now the integer variables; u variables mean on/off status
0765   if UC
0766     om.init_indexed_name('var', 'u', {nt});
0767     om.init_indexed_name('var', 'v', {nt});
0768     om.init_indexed_name('var', 'w', {nt});
0769     vt0 = char('B' * ones(1, ng));  % default variable type for u is binary
0770     for t = 1:nt
0771       umin = zeros(ng, 1);
0772       umax = ones(ng, 1);
0773       % min up/down restrictions on u
0774       if ~mdi.UC.CyclicCommitment
0775         % min up time has not passed yet since startup occured, force ON
0776         umin( (mdi.UC.InitialState > 0) & ...
0777               (t+mdi.UC.InitialState-mdi.UC.MinUp <= 0) ) = 1;
0778         % min down time has not passed yet since shutdown occured, force OFF
0779         umax( (mdi.UC.InitialState < 0) & ...
0780               (t-mdi.UC.InitialState-mdi.UC.MinDown <= 0) ) = 0;
0781       end
0782       % set limits for units forced ON or forced OFF
0783       iON = find(mdi.UC.CommitKey(:,t) == 2);
0784       iOFF = find(mdi.UC.CommitKey(:,t) == -1);
0785       umin(iON)  = 1;
0786       umax(iOFF) = 0;
0787 
0788       % set variable types
0789       vt = vt0;                 % initialize all variable types to binary
0790       vt(umin == umax) = 'C';   % make continuous for those that are fixed
0791       
0792       om.add_var('u', {t}, ng, zeros(ng, 1), umin, umax, vt);
0793     end
0794     % v variables mean startup events
0795     for t = 1:nt
0796       om.add_var('v', {t}, ng, zeros(ng, 1), zeros(ng, 1), ones(ng, 1));
0797     end
0798     % w variables mean shutdown events
0799     for t = 1:nt
0800       om.add_var('w', {t}, ng, zeros(ng, 1), zeros(ng, 1), ones(ng, 1));
0801     end
0802   end
0803   % An external program may be using coordination with AC flows, and in
0804   % that case we need corresponding Qg variables, whose only function is to
0805   % be constrained to zero if the commitment decision asks for that.
0806   if mdi.QCoordination
0807     om.init_indexed_name('var', 'Qg', {nt, nj_max, nc_max+1});
0808     for t = 1:nt
0809       for j = 1:mdi.idx.nj(t)
0810         for k = 1:mdi.idx.nc(t,j)+1
0811           mpc = mdi.flow(t,j,k).mpc;
0812           genmask = mpc.gen(:,GEN_STATUS) > 0;
0813           q0 = genmask .* mpc.gen(:, QG) / baseMVA;
0814           if UC     % relax bounds here, enforced by uQmax, uQmin constraints
0815             qmin = genmask .* (min(mpc.gen(:, QMIN) / baseMVA, 0) - 1);
0816             qmax = genmask .* (max(mpc.gen(:, QMAX) / baseMVA, 0) + 1);
0817           else      % enforce bounds here, subject to flow's GEN_STATUS
0818             qmin = genmask .* mpc.gen(:, QMIN) / baseMVA;
0819             qmax = genmask .* mpc.gen(:, QMAX) / baseMVA;
0820           end
0821           om.add_var('Qg', {t,j,k}, ng, q0, qmin, qmax);
0822         end
0823       end
0824     end
0825   end
0826   nvars = om.getN('var');
0827   mdi.idx.nvars = nvars;
0828 
0829   % Construct mechanism to keep track of expected storage states. This is
0830   % used for building some constraints in some types of problems, most
0831   % notably when there is a constraint on the terminal expected storage,
0832   % but it is also needed when there is a value associated with leftover
0833   % storage. Unfortunately this requires the construction of a substantial
0834   % mechanism for computing the expected terminal storage at any end-point
0835   % of the transition tree, be it a contingency or the terminal of the central
0836   % path at the end of the horizon.  Let SF(t) be the terminal storage value
0837   % at the end of the t-th period, assuming that we make it there without
0838   % contingencies, and let SI(t) be the initial amount of storage at that
0839   % same period. Then
0840   %     SI(t) = D(t) * SF(t-1).
0841   %     SF(t) = B1(t) * SI(t) + B2*[G(t)*x + H(t)]*x, and
0842   % Here, D(t) is created from the probability transition matrix, restricted
0843   % to transitions from base cases at t-1 to base cases at t. This allows
0844   % us to write a recursion. If the general form for SI(t),SF(t) is
0845   %     SI(t) = Li(t)*S0 + Mg(t)*x + Mh(t)*x,
0846   %     SF(t) = Lf(t)*S0 + Ng(t)*x + Nh(t)*x,
0847   % it turns out that the recursion is
0848   %     L(1) = D(1); Mg(1) = Mh(1) = 0; Ng(1) = G(1); Nh(1) = H(1);
0849   %     for t=2:nt
0850   %         Li(t) = D(t)*Lf(t-1) = D(t)*B1(t-1)*Li(t-1);
0851   %         Lf(t) = B1(t)*Li(t) = B1(t)*D(t)*Lf(t-1);
0852   %         Mg(t) = D(t)*Ng(t-1);
0853   %         Mh(t) = D(t)*Nh(t-1);
0854   %         Ng(t) = B1(t)*Mg(t) + B2(t)*G(t);
0855   %         Nh(t) = B1(t)*Mh(t) + B2(t)*H(t);
0856   %     end
0857   %
0858   % If SI,SF are organized first by blocks for each storage unit and within
0859   % the blocks by scenario, then the D matrix is simply made up by
0860   % repeating the str.tstep(t).TransMat matrix ns times in the diagonal
0861   % and then the columns weighted by the probabilities of the basecases at
0862   % t-1 given that we remained in basecases, and the rows are weighted by
0863   % the inverse of the probabilites of the scenarios at the beginning of t.
0864   % D(1) is special though, where each block is an nj(1) x 1 vector of ones.
0865   % The B matrices are formed by stacking appropriately sized diagonal
0866   % matrices for each storage unit along the diagonal, where each component
0867   % is simply the i-th element of beta times an identity matrix.
0868   if verbose
0869     fprintf('- Building expected storage-tracking mechanism.\n');
0870   end
0871   if ns
0872     % The following code assumes that no more variables will be added
0873     vv = om.get_idx();
0874     for t = 1:nt
0875       nsxnjt = ns*mdi.idx.nj(t);
0876       % Form G(t), H(t), B1(t), B2(t)
0877       G = sparse(nsxnjt, nvars);
0878       H = sparse(nsxnjt, nvars);
0879       B1 = sparse(nsxnjt, nsxnjt);
0880       B2 = sparse(nsxnjt, nsxnjt);
0881       for j = 1:mdi.idx.nj(t)
0882         ii  = ((1:ns)'-1)*mdi.idx.nj(t)+j;
0883         jj1 = (vv.i1.Psc(t,j,1):vv.iN.Psc(t,j,1))';
0884         jj2 = (vv.i1.Psd(t,j,1):vv.iN.Psd(t,j,1))';
0885         G = G + sparse(ii, jj1, -mdi.Delta_T  *  InEff(:,t), nsxnjt, nvars);
0886         H = H + sparse(ii, jj2, -mdi.Delta_T ./ OutEff(:,t), nsxnjt, nvars);
0887         B1 = B1 + sparse(ii, ii, beta1(:,t), nsxnjt, nsxnjt);
0888         B2 = B2 + sparse(ii, ii, beta2(:,t), nsxnjt, nsxnjt);
0889       end
0890       if t == 1
0891         % form Li, Lf, Mg, Mh, Ng, Nh, B1, B2 for t == 1
0892         jlist = [];
0893         for i=1:ns
0894           jlist = [ jlist; i*ones(mdi.idx.nj(t),1) ];
0895         end
0896         mdi.tstep(t).Li  = sparse((1:nsxnjt)', jlist, 1, nsxnjt, ns);
0897         mdi.tstep(t).Lf  = B1 * mdi.tstep(t).Li;
0898         mdi.tstep(t).Mg  = sparse(nsxnjt, nvars);   % Initial one is all zeros
0899         mdi.tstep(t).Mh  = sparse(nsxnjt, nvars);   % Initial one is all zeros
0900         mdi.tstep(t).Ng  = B2 * G;
0901         mdi.tstep(t).Nh  = B2 * H;
0902       else
0903         % Form D(t)
0904         D = sparse(nsxnjt, ns*mdi.idx.nj(t-1));
0905         p1 = mdi.CostWeights(1,1:mdi.idx.nj(t-1),t-1)';
0906         p1 = p1 / sum(p1);      % sigma(t)
0907         p2 = mdi.tstep(t).TransMat * p1;
0908         Di = spdiags(1./p2, 0, mdi.idx.nj(t), mdi.idx.nj(t)) * ...
0909                 sparse(mdi.tstep(t).TransMat) * ...
0910                 spdiags(p1, 0, mdi.idx.nj(t-1), mdi.idx.nj(t-1));
0911         for i = 1:ns
0912           D((i-1)*mdi.idx.nj(t)+1:i*mdi.idx.nj(t), (i-1)*mdi.idx.nj(t-1)+1:i*mdi.idx.nj(t-1)) = Di;
0913         end
0914         % Apply recursion, form Li, Lf, Mg, Mh, Ng, Nh
0915         mdi.tstep(t).Li = D  * mdi.tstep(t-1).Lf;
0916         mdi.tstep(t).Lf = B1 * mdi.tstep(t).Li;
0917         mdi.tstep(t).Mg = D * mdi.tstep(t-1).Ng;
0918         mdi.tstep(t).Mh = D * mdi.tstep(t-1).Nh;
0919         mdi.tstep(t).Ng = B1 * mdi.tstep(t).Mg + B2 * G;
0920         mdi.tstep(t).Nh = B1 * mdi.tstep(t).Mh + B2 * H;
0921       end
0922       mdi.tstep(t).G = G;
0923       mdi.tstep(t).H = H;
0924     end
0925   end
0926 
0927   % Now for the constraint indexing and creation.
0928   if verbose
0929     fprintf('- Building constraint submatrices.\n');
0930   end
0931   baseMVA = mdi.mpc.baseMVA;
0932   om.init_indexed_name('lin', 'Pmis', {nt, nj_max, nc_max+1});
0933   if mdi.DCMODEL
0934     % Construct all load flow equations using a DC flow model
0935     if verbose
0936       fprintf('  - Building DC flow constraints.\n');
0937     end
0938     om.init_indexed_name('lin', 'Pf', {nt, nj_max, nc_max+1});
0939     for t = 1:nt
0940       for j = 1:mdi.idx.nj(t)
0941         for k = 1:mdi.idx.nc(t,j)+1
0942           % First the flow constraints
0943           mpc = mdi.flow(t,j,k).mpc;
0944           ion = find(mpc.branch(:, BR_STATUS));
0945           [Bdc, Bl, Psh, PLsh] = makeBdc(baseMVA, mpc.bus, mpc.branch(ion,:));
0946           mdi.flow(t,j,k).PLsh = PLsh;     %% save for computing flows later
0947           negCg = sparse(mpc.gen(:,GEN_BUS), (1:ng)', -1, ...
0948                         mdi.idx.nb(t,j,k), ng);
0949           A = [Bdc negCg];
0950           b = -(mpc.bus(:,PD)+mpc.bus(:,GS))/baseMVA-Psh;
0951           vs = struct('name', {'Va', 'Pg'}, 'idx', {{t,j,k}, {t,j,k}});
0952           om.add_lin_constraint('Pmis', {t,j,k}, A, b, b, vs);
0953           % Then the thermal limits
0954           tmp = mpc.branch(ion,RATE_A)/baseMVA;
0955           iuncon = find(~tmp);
0956           tmp(iuncon) = Inf(size(iuncon));
0957           vs = struct('name', {'Va'}, 'idx', {{t,j,k}});
0958           om.add_lin_constraint('Pf', {t,j,k}, Bl, -tmp-PLsh, tmp-PLsh, vs);
0959         end
0960       end
0961     end
0962   else
0963     if verbose
0964       fprintf('  - Building load balance constraints.\n');
0965     end
0966     % Set simple generation - demand = 0 equations, one for each flow
0967     for t = 1:nt
0968       for j = 1:mdi.idx.nj(t)
0969         for k = 1:mdi.idx.nc(t,j)+1
0970           mpc = mdi.flow(t,j,k).mpc;
0971           A = sparse(ones(1, ng));
0972           b = 1.0*sum(mpc.bus(:, PD)+mpc.bus(:,GS))/baseMVA;
0973           vs = struct('name', {'Pg'}, 'idx', {{t,j,k}});
0974           om.add_lin_constraint('Pmis', {t,j,k}, A, b, b, vs);
0975         end
0976       end
0977     end
0978   end
0979   if mdi.IncludeFixedReserves
0980     if verbose
0981       fprintf('  - Building fixed zonal reserve constraints.\n');
0982     end
0983     om.init_indexed_name('lin', 'Pg_plus_R', {nt, nj_max, nc_max+1});
0984     om.init_indexed_name('lin', 'Rreq', {nt, nj_max, nc_max+1});
0985     for t = 1:nt
0986       for j = 1:mdi.idx.nj(t)
0987         for k = 1:mdi.idx.nc(t,j)+1
0988           % First the flow constraints
0989           mpc = mdi.flow(t,j,k).mpc;
0990           r = mdi.FixedReserves(t,j,k);
0991           ngr = length(r.igr);
0992           I = speye(ngr);
0993           Ar = sparse(1:ngr, r.igr, 1, ngr, ng);
0994           if UC
0995             A = [Ar I  ...
0996                   sparse(1:ngr, r.igr, -mpc.gen(r.igr, PMAX) / baseMVA, ngr, ng)];
0997             u = zeros(ngr, 1);
0998             vs = struct('name', {'Pg', 'R', 'u'}, 'idx', {{t,j,k}, {t,j,k}, {t}});
0999           else
1000             A = [Ar I];
1001             u = mpc.gen(r.igr, PMAX) / baseMVA;
1002             vs = struct('name', {'Pg', 'R'}, 'idx', {{t,j,k}, {t,j,k}});
1003           end
1004           om.add_lin_constraint('Pg_plus_R', {t,j,k}, A, [], u, vs);
1005           A = r.zones(:, r.igr);
1006           l = r.req / mpc.baseMVA;
1007           vs = struct('name', {'R'}, 'idx', {{t,j,k}});
1008           om.add_lin_constraint('Rreq', {t,j,k}, A, l, [], vs);
1009         end
1010       end
1011     end
1012   end
1013   
1014   % Set relationships between generator injections and charge/discharge
1015   % variables (-pg + psc + psd = 0)
1016   if verbose && ~isempty(mdi.Storage.UnitIdx)
1017     fprintf('  - Splitting storage injections into charge/discharge.\n');
1018   end
1019   om.init_indexed_name('lin', 'Ps', {nt, nj_max, nc_max+1});
1020   for t = 1:nt
1021     for j = 1:mdi.idx.nj(t)
1022       for k = 1:mdi.idx.nc(t,j)+1
1023         A = [sparse((1:ns)', mdi.Storage.UnitIdx, -1, ns, ng) Ins Ins];
1024         b = zeros(ns, 1);
1025         vs = struct('name', {'Pg', 'Psc', 'Psd'}, 'idx', {{t,j,k}, {t,j,k}, {t,j,k}});
1026         om.add_lin_constraint('Ps', {t,j,k}, A, b, b, vs);
1027       end
1028     end
1029   end
1030   
1031   % Construct y-variable restrictions on piecewise-linear costs. Note that
1032   % the restriction lines are computed using the full non-scaled cost in
1033   % gencost; any weighting of the cost must be then specified later in the
1034   % cost coefficients hitting the y variables (not 1 anymore).  Do it for
1035   % a complete c3sopf cell taking advantage of the fact that all p
1036   % injections are contiguous, as are all y variables for a c3sopf cell.
1037   % Also, note that makeAy assumes that every gen is online, which is
1038   % consistent with our formulation for unit commitment
1039   if verbose
1040     fprintf('  - Building CCV constraints for piecewise-linear costs.\n');
1041   end
1042   om.init_indexed_name('lin', 'ycon', {nt, nj_max, nc_max+1});
1043   for t = 1:nt,
1044     for j = 1:mdi.idx.nj(t)
1045       for k = 1:mdi.idx.nc(t,j)+1
1046         mpc = mdi.flow(t,j,k).mpc;
1047         [A, u] = makeAy(baseMVA, ng, mpc.gencost, 1, [], ng+1);
1048         vs = struct('name', {'Pg', 'y'}, 'idx', {{t,j,k}, {t,j,k}});
1049         om.add_lin_constraint('ycon', {t,j,k}, A, [], u, vs);
1050       end
1051     end
1052   end
1053 
1054   % The actual deviations from base flow must not exceed physical ramp rates
1055   % we'll get negative multiplier for right bound, fix when picking up
1056   % lambdas.
1057   % At issue: generators ousted in a contingency clearly violate this
1058   % transition; do not include constraints for these or for generators
1059   % whose commitment key is -1; either of these two possibilities will
1060   % result in a GEN_STATUS of 0, so we use that as the indicator.
1061   if verbose
1062     fprintf('  - Building contingency reserve constraints.\n');
1063   end
1064   om.init_indexed_name('lin', 'rampcont', {nt, nj_max, nc_max+1});
1065   for t =1:nt
1066     for j = 1:mdi.idx.nj(t)
1067       for k = 2:mdi.idx.nc(t,j)+1
1068         ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
1069         ngtmp = length(ii);
1070         A = sparse((1:ngtmp)', ii, 1, ngtmp, ng);
1071         u = mdi.flow(t,j,k).mpc.gen(ii,RAMP_10)/baseMVA;
1072         vs = struct('name', {'Pg', 'Pg'}, 'idx', {{t,j,1}, {t,j,k}});
1073         om.add_lin_constraint('rampcont', {t,j,k}, [-A A], -u, u, vs);
1074       end
1075     end
1076   end
1077   % The usual alpha-controlled equality of P0 and Pc does not make
1078   % sense when having many scenarios and hence many P0's .  Ditch.
1079   %
1080   % Positive reserve variables are larger than all increment variables in
1081   % all scenarios and flows of a given time slice 0 <= rpp - dpp; these
1082   % are the ones that set the price of reserves. Include all units that are
1083   % potentially committed.
1084   om.init_indexed_name('lin', 'dPpRp', {nt, nj_max, nc_max+1});
1085   for t = 1:nt
1086     for j = 1:mdi.idx.nj(t);
1087       for k = 1:mdi.idx.nc(t,j)+1
1088         ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
1089         ngtmp = length(ii);
1090         A = sparse((1:ngtmp)', ii, 1, ngtmp, ng);
1091         l = zeros(ngtmp, 1);
1092         vs = struct('name', {'dPp', 'Rpp'}, 'idx', {{t,j,k}, {t}});
1093         om.add_lin_constraint('dPpRp', {t,j,k}, [-A A], l, [], vs);
1094       end
1095     end
1096   end
1097   % Negative reserve variables are larger than all decrement variables in
1098   % all scenarios and flows of a given time slice  0 <= rpm - dpm; these
1099   % are the ones that set the price of reserves. Include all units that are
1100   % potentially committed.
1101   om.init_indexed_name('lin', 'dPmRm', {nt, nj_max, nc_max+1});
1102   for t = 1:nt
1103     for j = 1:mdi.idx.nj(t);
1104       for k = 1:mdi.idx.nc(t,j)+1
1105         ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
1106         ngtmp = length(ii);
1107         A = sparse((1:ngtmp)', ii, 1, ngtmp, ng);
1108         l = zeros(ngtmp, 1);
1109         vs = struct('name', {'dPm', 'Rpm'}, 'idx', {{t,j,k}, {t}});
1110         om.add_lin_constraint('dPmRm', {t,j,k}, [-A A], l, [], vs);
1111       end
1112     end
1113   end
1114   % The difference between the injection and the contract
1115   % is equal to the inc minus the dec: Ptjk - Ptc = dPp - dPm
1116   % Include all units that are potentially committed.
1117   om.init_indexed_name('lin', 'dPdef', {nt, nj_max, nc_max+1});
1118   for t = 1:nt
1119     for j = 1:mdi.idx.nj(t);
1120       for k = 1:mdi.idx.nc(t,j)+1
1121         ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
1122         ngtmp = length(ii);
1123         A = sparse((1:ngtmp)', ii, 1, ngtmp, ng);
1124         b = zeros(ngtmp, 1);
1125         vs = struct('name', {'Pg', 'Pc', 'dPp', 'dPm'}, ...
1126                     'idx', {{t,j,k}, {t}, {t,j,k}, {t,j,k}});
1127         om.add_lin_constraint('dPdef', {t,j,k}, [A -A -A A], b, b, vs);
1128       end
1129     end
1130   end
1131   
1132   % Go on to load following ramping restrictions.  Note that these
1133   % restrictions apply even if there is a change in the commitment status
1134   % of the generator.
1135   %
1136   % First, bound upward ramping reserves from below by all base-case
1137   % ramping possibilities, 0 <= rrp(t) -p(t+1)(j2)0 + p(t)(j1)0.  A ramping
1138   % reserve is needed at time t to be able to change to the needed dispatch at
1139   % time t+1.  An initial ramping reserve (t=0) would be data, not a
1140   % variable, and while ramping transitions from t=0 to t=1 should be
1141   % enforced, we do not allocate the reserve for t = 0.
1142   % Note: in the event that some future reserves are already locked in, we may
1143   %       not want to start at t = 1
1144   if verbose
1145     fprintf('  - Building ramping transitions and reserve constraints.\n');
1146   end
1147   om.init_indexed_name('lin', 'Rrp', {nt, nj_max, nj_max});
1148   % First, do from t=1:nt-1, since the last one is different and may not
1149   % even exist depending on the type of horizon
1150   for t = 1:nt-1
1151     for j1 = 1:mdi.idx.nj(t)      % j1 is at time t
1152       for j2 = 1:mdi.idx.nj(t+1)  % j2 is at time t+1
1153         if mdi.tstep(t+1).TransMask(j2,j1)
1154           A = [Ing -Ing Ing];
1155           l = zeros(ng, 1);
1156           vs = struct('name', {'Pg', 'Pg', 'Rrp'}, ...
1157                       'idx', {{t,j1,1}, {t+1,j2,1}, {t}});
1158           om.add_lin_constraint('Rrp', {t,j1,j2}, A, l, [], vs);
1159         end
1160       end
1161     end
1162   end
1163   % Now, pay special attention to a possible last type of ramping
1164   % constraint. If the horizon involves a terminal value at t=nt+1, then
1165   % this must also be enforced; in this case, additional ramping
1166   % reserves procured for t=nt must be defined.  If this
1167   % condition does not apply, then these reserves are not needed.
1168   if ~mdi.OpenEnded
1169     % pterminal <= rrp(nt) + p(nt,j1,0)
1170     for j1 = 1:mdi.idx.nj(nt)
1171       A = [Ing Ing];
1172       l = mdi.TerminalPg/baseMVA;
1173       vs = struct('name', {'Pg', 'Rrp'}, ...
1174                   'idx', {{nt,j1,1}, {nt}});
1175       om.add_lin_constraint('Rrp', {nt,j1,1}, A, l, [], vs);
1176     end
1177   end
1178   % Now on to downward ramping reserves.
1179   % First, bound downward ramping reserves from below by all base-case
1180   % ramping possibilities, 0 <= rrm(t) + p(t+1)j20 - p(t)j10
1181   om.init_indexed_name('lin', 'Rrm', {nt, nj_max, nj_max});
1182   % First, do from t=1:nt-1, since the last one is different and may not
1183   % even exist depending on the type of horizon
1184   for t = 1:nt-1
1185     for j1 = 1:mdi.idx.nj(t)      % j1 is at time t
1186       for j2 = 1:mdi.idx.nj(t+1)  % j2 is at time t+1
1187         if mdi.tstep(t+1).TransMask(j2,j1)
1188           A = [-Ing Ing Ing];
1189           l = zeros(ng, 1);
1190           vs = struct('name', {'Pg', 'Pg', 'Rrm'}, ...
1191                       'idx', {{t,j1,1}, {t+1,j2,1}, {t}});
1192           om.add_lin_constraint('Rrm', {t,j1,j2}, A, l, [], vs);
1193         end
1194       end
1195     end
1196   end
1197   % Now, pay special attention to a possible last type of ramping
1198   % constraint. If the horizon involves a terminal value at t=nt+1, then
1199   % this must also be enforced; in this case, additional ramping
1200   % reserves procured for t=nt must be defined.  If this
1201   % condition does not apply, then these reserves are not needed.
1202   if ~mdi.OpenEnded
1203     % -pterminal <= rrm(nt) - p(nt,j1,0)
1204     for j1 = 1:mdi.idx.nj(nt)
1205       A = [-Ing Ing];
1206       l = -mdi.TerminalPg/baseMVA;
1207       vs = struct('name', {'Pg', 'Rrm'}, ...
1208                   'idx', {{nt,j1,1}, {nt}});
1209       om.add_lin_constraint('Rrm', {nt,j1,1}, A, l, [], vs);
1210     end
1211   end
1212   %
1213   %
1214   % Now for the storage restrictions.
1215   if ns
1216     if verbose
1217       fprintf('  - Building storage constraints.\n');
1218     end
1219     % First bound sm(t) based on sm(t-1), with sm(1) being bound by the initial
1220     % data; this is for base case trajectories only
1221     om.init_indexed_name('lin', 'Sm', {nt, nj_max});
1222     if mdi.Storage.ForceCyclicStorage
1223       % sm(1) - beta1*s0 + beta2*Delta_T*[eta_c*psc(1,j,0) + (1/eta_d)*psd(1,j,0)] <= 0
1224       for j = 1:mdi.idx.nj(1)
1225         A = [ diagBeta2EtaIn1 diagBeta2overEtaOut1 Ins -spdiags(beta1(:,1), 0, ns, ns)];
1226         u = zeros(ns, 1);
1227         vs = struct('name', {'Psc', 'Psd', 'Sm', 'S0'}, 'idx', {{1,j,1}, {1,j,1}, {1}, {}});
1228         om.add_lin_constraint('Sm', {1,j}, A, [], u, vs);
1229       end
1230     else
1231       % sm(1) + beta2*Delta_T*[eta_c*psc(1,j,0) + (1/eta_d)*psd(1,j,0)] <= beta1*Initial/baseMVA
1232       for j = 1:mdi.idx.nj(1)
1233         A = [ diagBeta2EtaIn1 diagBeta2overEtaOut1 Ins ];
1234         u = beta1(:,1).*mdi.Storage.InitialStorageLowerBound/baseMVA;
1235         vs = struct('name', {'Psc', 'Psd', 'Sm'}, 'idx', {{1,j,1}, {1,j,1}, {1}});
1236         om.add_lin_constraint('Sm', {1,j}, A, [], u, vs);
1237       end
1238     end
1239     % Then the rest of the periods
1240     % sm(t) - beta1*(rho(t)*sm(t-1) + (1-rho(t))*s_I(t,j)) + beta2*Delta_T*[eta_c*psc(t,j,0) + (1/eta_d)*psd(t,j,0)] <= 0
1241     % where s_I(t,j) = L_I(t,j) * s0 + (Mg(t,j)+Mh(t,j)) * x
1242     for t = 2:nt
1243       for j = 1:mdi.idx.nj(t)
1244         Mj = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :) + ...
1245              mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1246         Lij = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1247         diag1minusRhoBeta1 = spdiags((1-rho(:,t)) .* beta1(:,t), 0, ns, ns);
1248         A = sparse([1:ns,1:ns,1:ns,1:ns]', ...
1249                    [vv.i1.Psc(t,j,1):vv.iN.Psc(t,j,1), vv.i1.Psd(t,j,1):vv.iN.Psd(t,j,1), vv.i1.Sm(t-1):vv.iN.Sm(t-1), vv.i1.Sm(t):vv.iN.Sm(t)]', ...
1250                    [beta2EtaIn(:,t); beta2overEtaOut(:,t); -beta1(:,t).*rho(:,t); ones(ns,1)], ...
1251                    ns, nvars) ...
1252                 - diag1minusRhoBeta1 * Mj;
1253         if mdi.Storage.ForceCyclicStorage
1254           As0 = sparse(ns, nvars);
1255           As0(:, vv.i1.S0:vv.iN.S0) = -diag1minusRhoBeta1 * Lij;
1256           A = A + As0;
1257           u = zeros(ns, 1);
1258         else
1259           u = full(diag1minusRhoBeta1 * Lij * mdi.Storage.InitialStorage/baseMVA);
1260         end
1261         om.add_lin_constraint('Sm', {t,j}, A, [], u);
1262       end
1263     end
1264     % Do the same we did for sm(t) for sp(t). First the initial step ...
1265     om.init_indexed_name('lin', 'Sp', {nt, nj_max});
1266     if mdi.Storage.ForceCyclicStorage
1267       % -sp(1) + beta1*s0 - beta2*Delta_T*[eta_c*psc(1,j,0) + (1/eta_d)*psd(1,j,0)] <= 0
1268       for j = 1:mdi.idx.nj(1)
1269         A = [ -diagBeta2EtaIn1 -diagBeta2overEtaOut1 -Ins spdiags(beta1(:,1), 0, ns, ns) ];
1270         u = zeros(ns, 1);
1271         vs = struct('name', {'Psc', 'Psd', 'Sp', 'S0'}, 'idx', {{1,j,1}, {1,j,1}, {1}, {}});
1272         om.add_lin_constraint('Sp', {1,j}, A, [], u, vs);
1273       end
1274     else
1275       % -sp(1) - beta2*Delta_T*[eta_c*psc(1,j,0) + (1/eta_d)*psd(1,j,0)] <= -beta1*Initial/baseMVA
1276       for j = 1:mdi.idx.nj(1)
1277         A = [ -diagBeta2EtaIn1 -diagBeta2overEtaOut1 -Ins ];
1278         u = -beta1(:,1).*mdi.Storage.InitialStorageUpperBound/baseMVA;
1279         vs = struct('name', {'Psc', 'Psd', 'Sp'}, 'idx', {{1,j,1}, {1,j,1}, {1}});
1280         om.add_lin_constraint('Sp', {1,j}, A, [], u, vs);
1281       end
1282     end
1283     % Then the rest of the periods
1284     % -sp(t) + beta1*(rho(t)*sp(t-1) + (1-rho(t))*s_I(t,j)) - beta2*Delta_T*[eta_c*psc(t,j,0) + (1/eta_d)*psd(t,j,0)] <= 0
1285     % where s_I(t,j) = L_I(t,j) * s0 + (Mg(t,j)+Mh(t,j)) * x
1286     for t = 2:nt
1287       for j = 1:mdi.idx.nj(t)
1288         Mj = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :) + ...
1289              mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1290         Lij = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1291         diag1minusRhoBeta1 = spdiags((1-rho(:,t)) .* beta1(:,t), 0, ns, ns);
1292         A = sparse([1:ns,1:ns,1:ns,1:ns]', ...
1293                    [vv.i1.Psc(t,j,1):vv.iN.Psc(t,j,1), vv.i1.Psd(t,j,1):vv.iN.Psd(t,j,1), vv.i1.Sp(t-1):vv.iN.Sp(t-1), vv.i1.Sp(t):vv.iN.Sp(t)]', ...
1294                    [-beta2EtaIn(:,t); -beta2overEtaOut(:,t); beta1(:,t).*rho(:,t); -ones(ns,1)], ...
1295                    ns, nvars) ...
1296                 + diag1minusRhoBeta1 * Mj;
1297         if mdi.Storage.ForceCyclicStorage
1298           As0 = sparse(ns, nvars);
1299           As0(:, vv.i1.S0:vv.iN.S0) = diag1minusRhoBeta1 * Lij;
1300           A = A + As0;
1301           u = zeros(ns, 1);
1302         else
1303           u = full(-diag1minusRhoBeta1 * Lij * mdi.Storage.InitialStorage/baseMVA);
1304         end
1305         om.add_lin_constraint('Sp', {t,j}, A, [], u);
1306       end
1307     end
1308     % Now go on and limit the amount of energy that can be used if a
1309     % contingency does happen. Bound sm first. First examine time period 1 wrt to initial
1310     % stored energy, then t=2 and on.
1311     om.init_indexed_name('lin', 'contSm', {nt, nj_max, nc_max+1});
1312     for j = 1:mdi.idx.nj(1)
1313       for k = 2:mdi.idx.nc(1,j)+1  %% NOTE NO k=1!!!
1314         if mdi.Storage.ForceCyclicStorage
1315           % beta4*Delta_T*[eta_c*psc(1,j,0) + (1/eta_d)*psd(1,j,0)] + beta3*Delta_T*[eta_c*psc(1,j,k) + (1/eta_d)*psd(1,j,k)] - beta5*s0 <= -sm_min(1)
1316           A = [ diagBeta4EtaIn1 diagBeta4overEtaOut1 diagBeta3EtaIn1 diagBeta3overEtaOut1 -spdiags(beta5(:,1), 0, ns, ns) ];
1317           u = -MinStorageLevel(:,1)/baseMVA;
1318           vs = struct('name', {'Psc', 'Psd', 'Psc', 'Psd', 'S0'}, ...
1319                       'idx', {{1,j,1}, {1,j,1}, {1,j,k}, {1,j,k}, {}});
1320         else
1321           % beta4*Delta_T*[eta_c*psc(1,j,0) + (1/eta_d)*psd(1,j,0)] + beta3*Delta_T*[eta_c*psc(1,j,k) + (1/eta_d)*psd(1,j,k)] <= beta5*Initial/baseMVA - sm_min(1)
1322           A = [ diagBeta4EtaIn1 diagBeta4overEtaOut1 diagBeta3EtaIn1 diagBeta3overEtaOut1 ];
1323           u = (beta5(:,1).*mdi.Storage.InitialStorageLowerBound - MinStorageLevel(:,1))/baseMVA;
1324           vs = struct('name', {'Psc', 'Psd', 'Psc', 'Psd'}, ...
1325                       'idx', {{1,j,1}, {1,j,1}, {1,j,k}, {1,j,k}});
1326         end
1327         om.add_lin_constraint('contSm', {1,j,k}, A, [], u, vs);
1328       end
1329     end
1330     % then the rest of the periods
1331     % -beta5*(rho(t)*sm(t-1) + (1-rho(t))*s_I(t,j)) + beta4*Delta_T*[eta_c*psc(t,j,0) + (1/eta_d)*psd(t,j,0)] + beta3*Delta_T*[eta_c*psc(t,j,k) + (1/eta_d)*psd(t,j,k)] <= -sm_min(t)
1332     % where s_I(t,j) = L_I(t,j) * s0 + (Mg(t,j)+Mh(t,j)) * x
1333     for t = 2:nt
1334       for j = 1:mdi.idx.nj(t)
1335         for k = 2:mdi.idx.nc(t,j)+1
1336           Mj = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :) + ...
1337                mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1338           Lij = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1339           diag1minusRhoBeta5 = spdiags((1-rho(:,t)) .* beta5(:,t), 0, ns, ns);
1340           A = sparse([1:ns,1:ns,1:ns,1:ns,1:ns]', ...
1341                      [vv.i1.Psc(t,j,1):vv.iN.Psc(t,j,1), vv.i1.Psd(t,j,1):vv.iN.Psd(t,j,1), vv.i1.Psc(t,j,k):vv.iN.Psc(t,j,k), vv.i1.Psd(t,j,k):vv.iN.Psd(t,j,k), vv.i1.Sm(t-1):vv.iN.Sm(t-1)]', ...
1342                      [beta4EtaIn(:,t); beta4overEtaOut(:,t); beta3EtaIn(:,t); beta3overEtaOut(:,t); -beta5(:,t).*rho(:,t)], ...
1343                      ns, nvars) ...
1344                   - diag1minusRhoBeta5 * Mj;
1345           u = -MinStorageLevel(:,t)/baseMVA;
1346           if mdi.Storage.ForceCyclicStorage
1347             As0 = sparse(ns, nvars);
1348             As0(:, vv.i1.S0:vv.iN.S0) = -diag1minusRhoBeta5 * Lij;
1349             A = A + As0;
1350           else
1351             u = u + diag1minusRhoBeta5 * Lij * mdi.Storage.InitialStorageLowerBound/baseMVA;
1352           end
1353           om.add_lin_constraint('contSm', {t,j,k}, A, [], u);
1354         end
1355       end
1356     end
1357     % Bound sp first. First examine time period 1 wrt to initial
1358     % stored energy, then t=2 and on.
1359     om.init_indexed_name('lin', 'contSp', {nt, nj_max, nc_max+1});
1360     for j = 1:mdi.idx.nj(1)
1361       for k = 2:mdi.idx.nc(1,j)+1
1362         if mdi.Storage.ForceCyclicStorage
1363           % -beta4*Delta_T*[eta_c*psc(1,j,0) + (1/eta_d)*psd(1,j,0)] - beta3*Delta_T*[eta_c*psc(1,j,k) + (1/eta_d)*psd(1,j,k)] + beta5*s0 <= sp_max(1)
1364           A = [ -diagBeta4EtaIn1 -diagBeta4overEtaOut1 -diagBeta3EtaIn1 -diagBeta3overEtaOut1 spdiags(beta5(:,1), 0, ns, ns)];
1365           u = MaxStorageLevel(:,1)/baseMVA;
1366           vs = struct('name', {'Psc', 'Psd', 'Psc', 'Psd', 'S0'}, ...
1367                       'idx', {{1,j,1}, {1,j,1}, {1,j,k}, {1,j,k}, {}});
1368         else
1369           % -beta4*Delta_T*[eta_c*psc(1,j,0) + (1/eta_d)*psd(1,j,0)] - beta3*Delta_T*[eta_c*psc(1,j,k) + (1/eta_d)*psd(1,j,k)] <= -beta5*Initial/baseMVA + sp_max(1)
1370           A = [ -diagBeta4EtaIn1 -diagBeta4overEtaOut1 -diagBeta3EtaIn1 -diagBeta3overEtaOut1 ];
1371           u = (MaxStorageLevel(:,1) - beta5(:,1).*mdi.Storage.InitialStorageUpperBound)/baseMVA;
1372           vs = struct('name', {'Psc', 'Psd', 'Psc', 'Psd'}, ...
1373                       'idx', {{1,j,1}, {1,j,1}, {1,j,k}, {1,j,k}});
1374         end
1375         om.add_lin_constraint('contSp', {1,j,k}, A, [], u, vs);
1376       end
1377     end
1378     % then the rest of the periods
1379     % beta5*(rho(t)*sp(t-1) + (1-rho(t))*s_I(t,j)) - beta4*Delta_T*[eta_c*psc(t,j,0) + (1/eta_d)*psd(t,j,0)] - beta3*Delta_T*[eta_c*psc(t,j,k) + (1/eta_d)*psd(t,j,k)] <= sp_max(t)
1380     % where s_I(t,j) = L_I(t,j) * s0 + (Mg(t,j)+Mh(t,j)) * x
1381     for t = 2:nt
1382       for j = 1:mdi.idx.nj(t)
1383         for k = 2:mdi.idx.nc(t,j)+1
1384           Mj = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :) + ...
1385                mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1386           Lij = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1387           diag1minusRhoBeta5 = spdiags((1-rho(:,t)) .* beta5(:,t), 0, ns, ns);
1388           A = sparse([1:ns,1:ns,1:ns,1:ns,1:ns]', ...
1389                      [vv.i1.Psc(t,j,1):vv.iN.Psc(t,j,1), vv.i1.Psd(t,j,1):vv.iN.Psd(t,j,1), vv.i1.Psc(t,j,k):vv.iN.Psc(t,j,k), vv.i1.Psd(t,j,k):vv.iN.Psd(t,j,k), vv.i1.Sp(t-1):vv.iN.Sp(t-1)]', ...
1390                      [-beta4EtaIn(:,t); -beta4overEtaOut(:,t); -beta3EtaIn(:,t); -beta3overEtaOut(:,t); beta5(:,t).*rho(:,t)], ...
1391                      ns, nvars) ...
1392                   + diag1minusRhoBeta5 * Mj;
1393           u = MaxStorageLevel(:,t)/baseMVA;
1394           if mdi.Storage.ForceCyclicStorage
1395             As0 = sparse(ns, nvars);
1396             As0(:, vv.i1.S0:vv.iN.S0) = diag1minusRhoBeta5 * Lij;
1397             A = A + As0;
1398           else
1399             u = u - diag1minusRhoBeta5 * Lij * mdi.Storage.InitialStorageUpperBound/baseMVA;
1400           end
1401           om.add_lin_constraint('contSp', {t,j,k}, A, [], u);
1402         end
1403       end
1404     end
1405   end
1406   
1407   % Now, if required, constrain the expected terminal storage quantity; two
1408   % different ways:
1409   if mdi.Storage.ForceExpectedTerminalStorage && mdi.Storage.ForceCyclicStorage
1410     error('most: ForceExpectedTerminalStorage and ForceCyclicStorage cannot be simultaneously true.');
1411   end
1412   if ns
1413     % The following code assumes that no more variables will be added
1414     if mdi.Storage.ForceExpectedTerminalStorage
1415       % 1) Constrain the expected terminal storage to be some target value
1416       A = sparse(ns, nvars);
1417       b = zeros(ns, 1);
1418       for j = 1:mdi.idx.nj(nt)
1419         Ngj = mdi.tstep(nt).Ng( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1420         Nhj = mdi.tstep(nt).Nh( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1421         Lfj = mdi.tstep(nt).Lf( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1422         A = A + mdi.CostWeights(1,j,nt) * (Ngj + Nhj);
1423         b = b + mdi.CostWeights(1,j,nt) * (Lfj * mdi.Storage.InitialStorage) / baseMVA;
1424       end
1425       endprob = sum(mdi.CostWeights(1,1:mdi.idx.nj(nt),nt)');
1426       A = (1/endprob) * A;
1427       b = (1/endprob) * b;
1428       l = mdi.Storage.ExpectedTerminalStorageMin / baseMVA - b;
1429       u = mdi.Storage.ExpectedTerminalStorageMax / baseMVA - b;
1430       om.add_lin_constraint('ESnt', A, l, u);
1431     elseif mdi.Storage.ForceCyclicStorage
1432       % 2) Constrain the initial storage (a variable) to be the same as the final expected storage
1433       A = sparse(ns, nvars);
1434       for j = 1:mdi.idx.nj(nt)
1435         Ngj = mdi.tstep(nt).Ng( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1436         Nhj = mdi.tstep(nt).Nh( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1437         A = A + mdi.CostWeights(1,j,nt) * (Ngj + Nhj);
1438       end
1439       endprob = sum(mdi.CostWeights(1,1:mdi.idx.nj(nt),nt)');
1440       A = (1/endprob) * A;
1441       for j = 1:mdi.idx.nj(nt)
1442         Lfj = mdi.tstep(nt).Lf(j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1443         A(:, vv.i1.S0:vv.iN.S0) = A(:, vv.i1.S0:vv.iN.S0) ...
1444                   + mdi.CostWeights(1,j,nt) * Lfj;
1445       end
1446       A(:, vv.i1.S0:vv.iN.S0) = (1/endprob) * A(:, vv.i1.S0:vv.iN.S0) - speye(ns);
1447       b = zeros(ns, 1);
1448       om.add_lin_constraint('ESnt', A, b, b);
1449     end
1450   end
1451 
1452   % Dynamical system contraints
1453   if nzds || nyds
1454     if verbose
1455       fprintf('  - Building dynamical system constraints.\n');
1456     end
1457     % Compute matrices that give the expected dispatch in time period t
1458     % given that we make it to that period, for all generators at once,
1459     % when multiplied by x, i.e. E[p(t)] = E(t) * x
1460     for t = 1:nt
1461       E = sparse(ng, nvars);
1462       for j = 1:mdi.idx.nj(t)
1463         for k = 1:mdi.idx.nc(t,j)+1;
1464           E = E + (mdi.CostWeightsAdj(k,j,t)/mdi.StepProb(t)) * sparse((1:ng)', (vv.i1.Pg(t,j,k):vv.iN.Pg(t,j,k))', 1, ng, nvars);
1465         end
1466       end
1467       mdi.tstep(t).E = E;
1468     end
1469   end
1470 
1471   % Form the dynamical system state equations and bound constraints on the
1472   % state vector
1473   if nzds
1474     om.init_indexed_name('lin', 'DSz', {ntds-1});
1475     b = zeros(nzds, 1);
1476     for t = 1:ntds-1
1477       if t <= nt  % We have p(t) available to drive the dynamical system up to t=nt
1478         % Form the constraint matrix, so B*E*x + A*z(t) - I*z(t+1) = 0
1479         A = mdi.dstep(t).B * mdi.tstep(t).E;
1480       else
1481         % The dynamical system horizon is longer than the injection planning
1482         % horizon and we don't know what p(t) is, but continue to drive the
1483         % dynamical system as if p(t) = 0 and perhaps take that into account
1484         % when setting Ymax, Ymin in this time window.  That is, A*z(t) - I*z(t+1) = 0
1485         A = sparse(nzds, nvars);
1486       end
1487       A(:, vv.i1.Z(t):vv.iN.Z(t)) = mdi.dstep(t).A;
1488       A(:, vv.i1.Z(t+1):vv.iN.Z(t+1)) = -speye(nzds);
1489       om.add_lin_constraint('DSz', {t}, A, b, b);
1490     end
1491   end
1492   
1493   % Form the output equations and their restrictions
1494   if nyds
1495     om.init_indexed_name('lin', 'DSy', {ntds});
1496     for t = 1:ntds
1497       if t <= nt
1498         A = mdi.dstep(t).D * mdi.tstep(t).E;
1499       else
1500         A = sparse(nyds, nvars);
1501       end
1502       if nzds
1503         A(:, vv.i1.Z(t):vv.iN.Z(t)) = mdi.dstep(t).C;
1504       end
1505       l = mdi.dstep(t).ymin;
1506       u = mdi.dstep(t).ymax;
1507       om.add_lin_constraint('DSy', {t}, A, l, u);
1508     end
1509   end
1510   
1511   % UNIT COMMITMENT
1512   if UC
1513     if verbose
1514       fprintf('  - Building unit commitment constraints.\n');
1515     end
1516     % u(t,i) - u(t-1,i) - v(t,i) + w(t,i) = 0
1517     om.init_indexed_name('lin', 'uvw', {nt});
1518     for t = 1:nt
1519       if t == 1
1520         % First for t=1 when u(t-1,i) is really u(0,i) or u(nt,i)
1521         if mdi.UC.CyclicCommitment
1522           vs = struct('name', {'u', 'u', 'v', 'w'}, 'idx', {{1}, {nt}, {1}, {1}});
1523           A = [Ing -Ing -Ing Ing];
1524           b = zeros(ng, 1);
1525         else
1526           vs = struct('name', {'u', 'v', 'w'}, 'idx', {{1}, {1}, {1}});
1527           A = [Ing -Ing Ing];
1528           b = (mdi.UC.InitialState > 0);
1529         end
1530       else
1531         % Then for rest of periods
1532         vs = struct('name', {'u', 'u', 'v', 'w'}, 'idx', {{t}, {t-1}, {t}, {t}});
1533         A = [Ing -Ing -Ing Ing];
1534         b = zeros(ng, 1);
1535       end
1536       om.add_lin_constraint('uvw', {t}, A, b, b, vs);
1537     end
1538     % Then continue with minimimum up time constraints. Again, two
1539     % different forms depending on whether the horizon is cyclical or not
1540     om.init_indexed_name('lin', 'minup', {nt, ng});
1541     for t = 1:nt
1542       for i = 1:ng
1543         ti = t-mdi.UC.MinUp(i)+1:t;
1544         if mdi.UC.CyclicCommitment     % window is circular
1545           for tt = 1:length(ti)
1546             if ti(tt) < 1
1547               ti(tt) = nt + ti(tt);
1548             end
1549           end
1550         end
1551         % limit to positive time
1552         % even with CyclicCommitment, in case MinUp is longer than horizon
1553         % (which implies always ON or always OFF)
1554         ti = ti(ti>0);
1555         vs = struct('name', {'u'}, 'idx', {{t}});
1556         A = sparse(1, i, -1, 1, ng);
1557         for tt = 1:length(ti)
1558             vs(end+1).name = 'v';
1559             vs(end).idx  = {ti(tt)};
1560             A = [A sparse(1, i, 1, 1, ng)];
1561         end
1562         om.add_lin_constraint('minup', {t, i}, A, [], 0, vs);
1563       end
1564     end
1565     % Continue with minimimum downtime constraints. Two
1566     % different forms depending on whether the horizon is cyclical or not
1567     om.init_indexed_name('lin', 'mindown', {nt, ng});
1568     for t = 1:nt
1569       for i = 1:ng
1570         ti = t-mdi.UC.MinDown(i)+1:t;
1571         if mdi.UC.CyclicCommitment     % window is circular
1572           for tt = 1:length(ti)
1573             if ti(tt) < 1
1574               ti(tt) = nt + ti(tt);
1575             end
1576           end
1577         end
1578         % limit to positive time
1579         % even with CyclicCommitment, in case MinDown is longer than horizon
1580         % (which implies always ON or always OFF)
1581         ti = ti(ti>0);
1582         vs = struct('name', {'u'}, 'idx', {{t}});
1583         A = sparse(1, i, 1, 1, ng);
1584         for tt = 1:length(ti)
1585             vs(end+1).name = 'w';
1586             vs(end).idx  = {ti(tt)};
1587             A = [A sparse(1, i, 1, 1, ng)];
1588         end
1589         om.add_lin_constraint('mindown', {t, i}, A, [], 1, vs);
1590       end
1591     end
1592     % Limit generation ranges based on commitment status; first Pmax;
1593     % p - u*Pmax <= 0
1594     % For contingent flows, however, if a generator is ousted as a result
1595     % of the contingency, then this constraint should not be enforced.
1596     om.init_indexed_name('lin', 'uPmax', {nt, nj_max, nc_max+1});
1597     for t = 1:nt
1598       for j = 1:mdi.idx.nj(t)
1599         for k = 1:mdi.idx.nc(t,j)+1
1600           mpc = mdi.flow(t,j,k).mpc;
1601           ii = find(mpc.gen(:, GEN_STATUS));
1602           nii = length(ii);
1603           vs = struct('name', {'Pg', 'u'}, 'idx', {{t,j,k}, {t}});
1604           A = [ sparse(1:nii, ii, 1, nii, ng) ...
1605                 sparse(1:nii, ii, -mpc.gen(ii, PMAX)/baseMVA, nii, ng) ];
1606           u = zeros(nii, 1);
1607           om.add_lin_constraint('uPmax', {t,j,k}, A, [], u, vs);
1608         end
1609       end
1610     end
1611     % Then Pmin,  -p + u*Pmin <= 0
1612     om.init_indexed_name('lin', 'uPmin', {nt, nj_max, nc_max+1});
1613     for t = 1:nt
1614       for j = 1:mdi.idx.nj(t)
1615         for k = 1:mdi.idx.nc(t,j)+1
1616           mpc = mdi.flow(t,j,k).mpc;
1617           ii = find(mpc.gen(:, GEN_STATUS));
1618           nii = length(ii);
1619           vs = struct('name', {'Pg', 'u'}, 'idx', {{t,j,k}, {t}});
1620           A = [ sparse(1:nii, ii, -1, nii, ng) ...
1621                 sparse(1:nii, ii, mpc.gen(ii, PMIN)/baseMVA, nii, ng) ];
1622           u = zeros(nii, 1);
1623           om.add_lin_constraint('uPmin', {t,j,k}, A, [], u, vs);
1624         end
1625       end
1626     end
1627     % Then, if there is Qg coordination, do the same for Qg
1628     % q - u*Qmax <= 0
1629     % For contingent flows, however, if a generator is ousted as a result
1630     % of the contingency, then this constraint should not be enforced.
1631     if mdi.QCoordination
1632       om.init_indexed_name('lin', 'uQmax', {nt, nj_max, nc_max+1});
1633       for t = 1:nt
1634         for j = 1:mdi.idx.nj(t)
1635           for k = 1:mdi.idx.nc(t,j)+1
1636             mpc = mdi.flow(t,j,k).mpc;
1637             ii = find(mpc.gen(:, GEN_STATUS));
1638             nii = length(ii);
1639             vs = struct('name', {'Qg', 'u'}, 'idx', {{t,j,k}, {t}});
1640             A = [ sparse(1:nii, ii, 1, nii, ng) ...
1641                   sparse(1:nii, ii, -mpc.gen(ii, QMAX)/baseMVA, nii, ng) ];
1642             u = zeros(nii, 1);
1643             om.add_lin_constraint('uQmax', {t,j,k}, A, [], u, vs);
1644           end
1645         end
1646       end
1647       % Then Qmin,  -q + u*Qmin <= 0
1648       om.init_indexed_name('lin', 'uQmin', {nt, nj_max, nc_max+1});
1649       for t = 1:nt
1650         for j = 1:mdi.idx.nj(t)
1651           for k = 1:mdi.idx.nc(t,j)+1
1652             mpc = mdi.flow(t,j,k).mpc;
1653             ii = find(mpc.gen(:, GEN_STATUS));
1654             nii = length(ii);
1655             vs = struct('name', {'Qg', 'u'}, 'idx', {{t,j,k}, {t}});
1656             A = [ sparse(1:nii, ii, -1, nii, ng) ...
1657                   sparse(1:nii, ii, mpc.gen(ii, QMIN)/baseMVA, nii, ng) ];
1658             u = zeros(nii, 1);
1659             om.add_lin_constraint('uQmin', {t,j,k}, A, [], u, vs);
1660           end
1661         end
1662       end
1663     end
1664   end
1665 
1666   if verbose
1667     fprintf('- Building cost structures.\n');
1668   end
1669   % Start building the cost.  Two main components, the input data cost and
1670   % the coordination cost are specified.  The coordination cost is assumed to
1671   % have been buit with knowledge of the variable structure, and is simply
1672   % passed on.  The input data cost is assembled into the appropriate
1673   % spots.
1674   %
1675   % f = 0.5 * x' * (H1 + Hcoord) * x + (C1' + Ccoord) * x + c1 + ccoord
1676 
1677   % First assign the ramping costs; H1 has few coefficients initially and
1678   % this should make the shuffling and reordering of coefficients more
1679   % efficient.  All other accesses to H1 will be diagonal insertions, which
1680   % take less time than anti-diagonal insertions.
1681   % First do first period wrt to InitialPg.
1682   if mdi.OpenEnded
1683     om.init_indexed_name('qdc', 'RampWear', {nt, nj_max, nj_max});
1684   else
1685     om.init_indexed_name('qdc', 'RampWear', {nt+1, nj_max, nj_max});
1686   end
1687   for j = 1:mdi.idx.nj(1)
1688     w = mdi.tstep(1).TransMat(j,1);  % the probability of going from initial state to jth
1689     Q = spdiags(w * baseMVA^2 * mdi.RampWearCostCoeff(:,1), 0, ng, ng);
1690     c = -w * baseMVA * mdi.RampWearCostCoeff(:,1) .* mdi.InitialPg;
1691     vs = struct('name', {'Pg'}, 'idx', {{1,j,1}});
1692     k0 = w * 0.5 * mdi.RampWearCostCoeff(:,1)' * mdi.InitialPg.^2;
1693     om.add_quad_cost('RampWear', {1,j,1}, Q, c, k0, vs);
1694   end
1695   % Then the remaining periods
1696   for t = 2:nt
1697     for j2 = 1:mdi.idx.nj(t)
1698       for j1 = 1:mdi.idx.nj(t-1)
1699         w = mdi.tstep(t).TransMat(j2,j1) * mdi.CostWeights(1, j1, t-1);
1700         h = w * baseMVA^2 * mdi.RampWearCostCoeff(:,t);
1701         i = (1:ng)';
1702         j = ng+(1:ng)';
1703         Q = sparse([i;j;i;j], [i;i;j;j], [h;-h;-h;h], 2*ng, 2*ng);
1704         vs = struct('name', {'Pg', 'Pg'}, 'idx', {{t-1,j1,1}, {t,j2,1}});
1705         om.add_quad_cost('RampWear', {t,j1,j2}, Q, zeros(2*ng,1), 0, vs);
1706       end
1707     end
1708   end
1709   % Finally, if there is a terminal state problem, apply cost to
1710   % the transition starting from t=nt.  Note that in this case
1711   % mdi.tstep(nt+1).TransMat must be defined! it is the only piece of data
1712   % that makes sense for nt+1; all other fields in mdi.tstep(nt+1) can be empty.
1713   if ~mdi.OpenEnded
1714     for j = 1:mdi.idx.nj(nt)
1715       w = mdi.tstep(nt+1).TransMat(1, j) * mdi.CostWeights(1, j, nt);
1716       Q = spdiags(w * baseMVA^2 * mdi.RampWearCostCoeff(:,nt+1), 0, ng, ng);
1717       c = -w * baseMVA * mdi.RampWearCostCoeff(:,nt+1) .* mdi.TerminalPg;
1718       vs = struct('name', {'Pg'}, 'idx', {{nt,j,1}});
1719       k0 = w * 0.5 * mdi.RampWearCostCoeff(:,nt+1)' * mdi.TerminalPg.^2;
1720       om.add_quad_cost('RampWear', {nt+1,j,1}, Q, c, k0, vs);
1721     end
1722   end
1723 
1724   % Now go on and assign energy, inc/dec and contingency reserves
1725   % costs for all committed units.
1726   om.init_indexed_name('qdc', 'Cp', {nt, nj_max, nc_max+1});
1727   om.init_indexed_name('qdc', 'Cy', {nt, nj_max, nc_max+1});
1728   om.init_indexed_name('qdc', 'Cpp', {nt, nj_max, nc_max+1});
1729   om.init_indexed_name('qdc', 'Cpm', {nt, nj_max, nc_max+1});
1730   if mdi.IncludeFixedReserves
1731     om.init_indexed_name('qdc', 'Rcost', {nt, nj_max, nc_max+1});
1732   end
1733   om.init_indexed_name('qdc', 'Crpp', {nt});
1734   om.init_indexed_name('qdc', 'Crpm', {nt});
1735   for t = 1:nt
1736     for j = 1:mdi.idx.nj(t)
1737       for k = 1:mdi.idx.nc(t,j)+1
1738         w = mdi.CostWeightsAdj(k,j,t);     %% NOTE (k,j,t) order !!!
1739 
1740         % weighted polynomial energy costs for committed units
1741         gc = mdi.flow(t,j,k).mpc.gencost;
1742         ipol = find(gc(:, MODEL) == POLYNOMIAL);
1743         if ~isempty(ipol)
1744           ncost = gc(ipol(1), NCOST);
1745           if all(gc(ipol, NCOST) == ncost)    %% uniform order of polynomials
1746             %% use vectorized code
1747             if ncost > 3
1748               error('most: polynomial generator costs of order higher than quadratic not supported');
1749             elseif ncost == 3
1750               Q = sparse(ipol, ipol, 2 * w * baseMVA^2*gc(ipol, COST), ng, ng);
1751             else
1752               Q = sparse(ng,ng);
1753             end
1754             c = zeros(ng, 1);
1755             if ncost >= 2
1756               c(ipol) = w * baseMVA*gc(ipol, COST+ncost-2);
1757             end
1758             k0 = w * sum(gc(ipol, COST+ncost-1));
1759           else                                %% non-uniform order of polynomials
1760             %% use a loop
1761             Q = sparse(ng,ng);
1762             c = zeros(ng, 1);
1763             for i = ipol'
1764               ncost = gc(i, NCOST);
1765               if ncost > 3
1766                 error('most: polynomial generator costs of order higher than quadratic not supported');
1767               elseif ncost == 3
1768                 Q(i,i) = 2 * w * baseMVA^2*gc(i, COST);
1769               end
1770               if ncost >= 2
1771                 c(i) = w * baseMVA*gc(i, COST+ncost-2);
1772               end
1773               k0 = w * gc(i, COST+ncost-1);
1774             end
1775           end
1776           vs = struct('name', {'Pg'}, 'idx', {{t,j,k}});
1777           om.add_quad_cost('Cp', {t,j,k}, Q, c, k0, vs);
1778         end
1779 
1780         % weighted y-variables for piecewise linear energy costs for committed units
1781         % ipwl = find( (mdi.flow(t,j,k).mpc.gen(:,GEN_STATUS) > 0) & (gc(:,MODEL) == PW_LINEAR));
1782         if mdi.idx.ny(t,j,k)
1783           c = w * ones(mdi.idx.ny(t,j,k),1);
1784           vs = struct('name', {'y'}, 'idx', {{t,j,k}});
1785           om.add_quad_cost('Cy', {t,j,k}, [], c, 0, vs);
1786         end
1787 
1788         % inc and dec offers for each flow
1789         c = w * baseMVA * mdi.offer(t).PositiveActiveDeltaPrice(:);
1790         vs = struct('name', {'dPp'}, 'idx', {{t,j,k}});
1791         om.add_quad_cost('Cpp', {t,j,k}, [], c, 0, vs);
1792         c = w * baseMVA * mdi.offer(t).NegativeActiveDeltaPrice(:);
1793         vs = struct('name', {'dPm'}, 'idx', {{t,j,k}});
1794         om.add_quad_cost('Cpm', {t,j,k}, [], c, 0, vs);
1795 
1796         % weighted fixed reserves cost
1797         if mdi.IncludeFixedReserves
1798           c = w * mdi.FixedReserves(t,j,k).cost(r.igr) * baseMVA;
1799           vs = struct('name', {'R'}, 'idx', {{t,j,k}});
1800           om.add_quad_cost('Rcost', {t,j,k}, [], c, 0, vs);
1801         end
1802       end
1803     end
1804     
1805     % contingency reserve costs
1806     c = baseMVA * mdi.StepProb(t) * mdi.offer(t).PositiveActiveReservePrice(:);
1807     vs = struct('name', {'Rpp'}, 'idx', {{t}});
1808     om.add_quad_cost('Crpp', {t}, [], c, 0, vs);
1809     c = baseMVA * mdi.StepProb(t) * mdi.offer(t).NegativeActiveReservePrice(:);
1810     vs = struct('name', {'Rpm'}, 'idx', {{t}});
1811     om.add_quad_cost('Crpm', {t}, [], c, 0, vs);
1812   end
1813   % Assign load following ramp reserve costs.  Do first nt-1 periods first
1814   om.init_indexed_name('qdc', 'Crrp', {mdi.idx.ntramp});
1815   om.init_indexed_name('qdc', 'Crrm', {mdi.idx.ntramp});
1816   for t = 1:nt-1,
1817     c = baseMVA * mdi.StepProb(t+1) * mdi.offer(t).PositiveLoadFollowReservePrice(:);
1818     vs = struct('name', {'Rrp'}, 'idx', {{t}});
1819     om.add_quad_cost('Crrp', {t}, [], c, 0, vs);
1820     c = baseMVA * mdi.StepProb(t+1) * mdi.offer(t).NegativeLoadFollowReservePrice(:);
1821     vs = struct('name', {'Rrm'}, 'idx', {{t}});
1822     om.add_quad_cost('Crrm', {t}, [], c, 0, vs);
1823   end
1824   % Then do last period if needed Terminal state case
1825   if ~mdi.OpenEnded
1826     %% are these costs missing a mdi.StepProb(t)?  -- rdz
1827     c = baseMVA * mdi.offer(nt).PositiveLoadFollowReservePrice(:);
1828     vs = struct('name', {'Rrp'}, 'idx', {{nt}});
1829     om.add_quad_cost('Crrp', {nt}, [], c, 0, vs);
1830     c = baseMVA * mdi.offer(nt).NegativeLoadFollowReservePrice(:);
1831     vs = struct('name', {'Rrm'}, 'idx', {{nt}});
1832     om.add_quad_cost('Crrm', {nt}, [], c, 0, vs);
1833   end
1834   % Assign startup/shutdown costs, if any, and fixed operating costs
1835   if UC
1836     om.init_indexed_name('qdc', 'c00', {nt});
1837     om.init_indexed_name('qdc', 'startup', {nt});
1838     om.init_indexed_name('qdc', 'shutdown', {nt});
1839     for t = 1:nt
1840       ww = zeros(ng, 1);
1841       for j = 1:mdi.idx.nj(t)
1842         for k = 1:mdi.idx.nc(t)+1
1843           ww = ww + mdi.CostWeightsAdj(k,j,t) * mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS);
1844         end
1845       end
1846       c = ww.*mdi.UC.c00(:,t);
1847       vs = struct('name', {'u'}, 'idx', {{t}});
1848       om.add_quad_cost('c00', {t}, [], c, 0, vs);
1849       c = mdi.StepProb(t)*mdi.flow(t,1,1).mpc.gencost(:, STARTUP);
1850       vs = struct('name', {'v'}, 'idx', {{t}});
1851       om.add_quad_cost('startup', {t}, [], c, 0, vs);
1852       c = mdi.StepProb(t)*mdi.flow(t,1,1).mpc.gencost(:, SHUTDOWN);
1853       vs = struct('name', {'w'}, 'idx', {{t}});
1854       om.add_quad_cost('shutdown', {t}, [], c, 0, vs);
1855     end
1856   end
1857   % Finally, assign any value to leftover stored energy
1858   if ns
1859     A1 = sparse(ns, ns);
1860     A2 = sparse(ns, nvars);
1861     A3 = sparse(ns, nvars);
1862     A4 = sparse(ns, nvars);
1863     A5 = sparse(ns, nvars);
1864     A6 = sparse(ns, nvars);
1865     A7 = sparse(ns, nvars);
1866     % The following code assumes that no more variables will be added
1867     vv = om.get_idx();
1868     for t = 1:nt
1869       % Compute cost coefficients for value of expected leftover storage
1870       % after a contingency
1871       for j = 1:mdi.idx.nj(t)
1872         % pick rows for jth base injections
1873         Gtj0  = mdi.tstep(t).G(  j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1874         Htj0  = mdi.tstep(t).H(  j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1875         Litj0 = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1876         Mgtj0 = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1877         Mhtj0 = mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1878         sum_psi_tjk = sum(mdi.CostWeights(2:mdi.idx.nc(t,j)+1,j,t));
1879         if t == nt
1880           A1 = A1 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta1(:,t), 0, ns, ns) * Litj0;
1881           A2 = A2 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta1(:,t), 0, ns, ns) * Mgtj0;
1882           A3 = A3 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta1(:,t), 0, ns, ns) * Mhtj0;
1883           A4 = A4 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta2(:,t), 0, ns, ns) * Gtj0;
1884           A5 = A5 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta2(:,t), 0, ns, ns) * Htj0;
1885         end
1886         A1 = A1 + sum_psi_tjk * spdiags(OutEff(:,t) .* beta5(:,t), 0, ns, ns) * Litj0;
1887         A2 = A2 + sum_psi_tjk * (spdiags(OutEff(:,t) .* beta5(:,t), 0, ns, ns) * Mgtj0 + spdiags(OutEff(:,t) .* beta4(:,t), 0, ns, ns) * Gtj0);
1888         A3 = A3 + sum_psi_tjk * (spdiags(OutEff(:,t) .* beta5(:,t), 0, ns, ns) * Mhtj0 + spdiags(OutEff(:,t) .* beta4(:,t), 0, ns, ns) * Htj0);
1889         for k = 2:mdi.idx.nc(t,j)+1
1890           ii  = (1:ns)';
1891           jj1 = (vv.i1.Psc(t,j,k):vv.iN.Psc(t,j,k))';
1892           jj2 = (vv.i1.Psd(t,j,k):vv.iN.Psd(t,j,k))';
1893           Gtjk = sparse(ii, jj1, -mdi.Delta_T  *  InEff(:,t), ns, nvars);
1894           Htjk = sparse(ii, jj2, -mdi.Delta_T ./ OutEff(:,t), ns, nvars);
1895           A6 = A6 + mdi.CostWeights(k,j,t) * spdiags(OutEff(:,t) .* beta3(:,t), 0, ns, ns) * Gtjk;
1896           A7 = A7 + mdi.CostWeights(k,j,t) * spdiags(OutEff(:,t) .* beta3(:,t), 0, ns, ns) * Htjk;
1897         end
1898       end
1899     end
1900     Cfstor = -baseMVA * ...
1901        (mdi.Storage.TerminalStoragePrice'      * (A2 + A3) + ...
1902         mdi.Storage.TerminalChargingPrice0'    * A4 + ...
1903         mdi.Storage.TerminalDischargingPrice0' * A5 + ...
1904         mdi.Storage.TerminalChargingPriceK'    * A6 + ...
1905         mdi.Storage.TerminalDischargingPriceK' * A7);
1906     if mdi.Storage.ForceCyclicStorage
1907       % If the horizon model for the storage is cyclic and therefore s0 is a
1908       % variable, then that initial storage must come at a cost,
1909       % (InitialStorageCost) otherwise the optimizer will force the gratis
1910       % s0 up just to have (possibly) more storage left at the end.
1911       Cfstor(vv.i1.S0:vv.iN.S0) = ...
1912         Cfstor(vv.i1.S0:vv.iN.S0) + ...
1913             baseMVA * mdi.Storage.InitialStorageCost';
1914       % and the term in the final expected storage related to s0 is also
1915       % not constant, so must be included in the objective function
1916       Cfstor(vv.i1.S0:vv.iN.S0) = ...
1917           Cfstor(vv.i1.S0:vv.iN.S0) - ...
1918           baseMVA * mdi.Storage.TerminalStoragePrice' * A1;
1919     end
1920     om.add_quad_cost('fstor', [], Cfstor', 0);
1921 
1922     % The following is a hack to make the storage state bounds tight;
1923     % assign them a very small cost
1924     om.init_indexed_name('qdc', 'SpSmFudge', {nt});
1925     c = 1e-2 * [-ones(ns,1); ones(ns,1)];
1926     for t = 1:nt
1927       vs = struct('name', {'Sm', 'Sp'}, 'idx', {{t}, {t}});
1928       om.add_quad_cost('SpSmFudge', {t}, [], c, 0, vs);
1929     end
1930   else
1931     Cfstor = sparse(1, nvars);
1932   end
1933 
1934   % Plug into struct
1935   if verbose
1936     fprintf('- Assembling full set of costs.\n');
1937   end
1938   [Q, c, k0] = om.params_quad_cost();
1939   mdi.QP.Cfstor = Cfstor;
1940   mdi.QP.H1 = Q;
1941   mdi.QP.C1 = c;
1942   mdi.QP.c1 = k0;
1943 end     % if mpopt.most.build_model
1944 
1945 % With all pieces of the cost in place, can proceed to build the total
1946 % cost now.
1947 mdi.QP.H = mdi.QP.H1;
1948 mdi.QP.C = mdi.QP.C1;
1949 mdi.QP.c = mdi.QP.c1;
1950 if isfield(mdi, 'CoordCost') && ...
1951         (~isempty(mdi.CoordCost.Cuser) || ~isempty(mdi.CoordCost.Huser))
1952   if verbose
1953     fprintf('- Adding coordination cost to standard cost.\n');
1954   end
1955   nvuser = length(mdi.CoordCost.Cuser);
1956   nvars = mdi.idx.nvars;
1957   mdi.QP.H = mdi.QP.H + ...
1958             [ mdi.CoordCost.Huser       sparse(nvuser,nvars-nvuser) ;
1959             sparse(nvars-nvuser,nvuser)  sparse(nvars-nvuser,nvars-nvuser) ];
1960   mdi.QP.C(1:nvuser) = mdi.QP.C(1:nvuser) +  mdi.CoordCost.Cuser(:);
1961   mdi.QP.c = mdi.QP.c + mdi.CoordCost.cuser;
1962   
1963 %   cp = struct('Cw', mdi.CoordCost.Cuser(:), ...
1964 %         'H', [ mdi.CoordCost.Huser     sparse(nvuser,nvars-nvuser) ;
1965 %             sparse(nvars-nvuser,nvuser) sparse(nvars-nvuser,nvars-nvuser) ]);
1966 %   om.add_legacy_cost('CoordCost', cp);
1967 end
1968 
1969 mdi.om = om;
1970 [vv, ll] = om.get_idx();
1971 if verbose
1972   fprintf('- Assembling full set of constraints.\n');
1973 end
1974 [mdi.QP.A, mdi.QP.l, mdi.QP.u] = om.params_lin_constraint();
1975 if verbose
1976   fprintf('- Assembling full set of variable bounds.\n');
1977 end
1978 [mdi.QP.x0, mdi.QP.xmin, mdi.QP.xmax, mdi.QP.vtype] = om.params_var();
1979 
1980 et_setup = toc(t0);
1981 t0 = tic;
1982 
1983 % Call solver!
1984 mdo = mdi;
1985 if mpopt.most.solve_model
1986   %% check consistency of model options (in case mdi was built in previous call)
1987   if mdi.DCMODEL ~= mo.DCMODEL
1988     error('MDI.DCMODEL inconsistent with MPOPT.most.dc_model');
1989   end
1990   if mdi.IncludeFixedReserves ~= mo.IncludeFixedReserves
1991     error('MDI.IncludeFixedReserves inconsistent with MPOPT.most.fixed_res (and possible presence of MDI.FixedReserves(t,j,k))');
1992   end
1993   if mdi.SecurityConstrained ~= mo.SecurityConstrained
1994     error('MDI.SecurityConstrained inconsistent with MPOPT.most.security_constraints (and possible presence of MDI.cont(t,j).contab)');
1995   end
1996   if mdi.QCoordination ~= mo.QCoordination
1997     error('MDI.QCoordination inconsistent with MPOPT.most.q_coordination');
1998   end
1999   if mdi.Storage.ForceCyclicStorage ~= mo.ForceCyclicStorage
2000     error('MDI.Storage.ForceCyclicStorage inconsistent with MPOPT.most.storage.cyclic');
2001   end
2002   if mdi.Storage.ForceExpectedTerminalStorage ~= mo.ForceExpectedTerminalStorage
2003     error('MDI.Storage.ForceExpectedTerminalStorage inconsistent with MPOPT.most.storage.terminal_target (and possible presence of MDI.Storage.ExpectedTerminalStorageAim|Min|Max)');
2004   end
2005   if mdi.UC.run ~= UC
2006     error('MDI.UC.run inconsistent with MPOPT.most.uc.run (and possible presence of MDI.UC.CommitKey)');
2007   end
2008   %% set options
2009   if any(any(mdi.QP.H))
2010     model = 'QP';
2011   else
2012     model = 'LP';
2013   end
2014   if UC
2015     model = ['MI' model];
2016   end
2017   mdo.QP.opt = mpopt2qpopt(mpopt, model, 'most');
2018   if verbose
2019     fprintf('- Calling %s solver.\n\n', model);
2020     fprintf('============================================================================\n\n');
2021   end
2022   if UC
2023     [mdo.QP.x, mdo.QP.f, mdo.QP.exitflag, mdo.QP.output, ...
2024             mdo.QP.lambda ] = miqps_matpower( mdi.QP.H, mdi.QP.C, ...
2025                 mdi.QP.A, mdi.QP.l, mdi.QP.u, mdi.QP.xmin, mdi.QP.xmax, ...
2026                 [], mdi.QP.vtype, mdo.QP.opt);
2027   else
2028     [mdo.QP.x, mdo.QP.f, mdo.QP.exitflag, mdo.QP.output, ...
2029             mdo.QP.lambda ] = qps_matpower( mdi.QP.H, mdi.QP.C, ...
2030                 mdi.QP.A, mdi.QP.l, mdi.QP.u, mdi.QP.xmin, mdi.QP.xmax, ...
2031                 [], mdo.QP.opt);
2032   end
2033   if mdo.QP.exitflag > 0
2034     success = 1;
2035     if verbose
2036       fprintf('\n============================================================================\n');
2037       fprintf('- MOST: %s solved successfully.\n', model);
2038     end
2039   else
2040     success = 0;
2041     if verbose
2042       fprintf('\n============================================================================\n');
2043       fprintf('- MOST: %s solver ''%s'' failed with exit flag = %d\n', model, mdo.QP.opt.alg, mdo.QP.exitflag);
2044     end
2045 %     fprintf('\n============================================================================\n');
2046 %     fprintf('- MOST: %s solver ''%s'' failed with exit flag = %d\n', model, mdo.QP.opt.alg, mdo.QP.exitflag);
2047 %     fprintf('  You can query the workspace to debug.\n')
2048 %     fprintf('  When finished, type the word "dbcont" to continue.\n\n');
2049 %     keyboard;
2050   end
2051   % Unpack results
2052   if verbose
2053     fprintf('- Post-processing results.\n');
2054   end
2055   if success
2056     for t = 1:nt
2057       if UC
2058         mdo.UC.CommitSched(:, t) = mdo.QP.x(vv.i1.u(t):vv.iN.u(t));
2059       end
2060       for j = 1:mdi.idx.nj(t)
2061         for k = 1:mdi.idx.nc(t,j)+1
2062           mpc = mdo.flow(t,j,k).mpc;      %% pull mpc from output struct
2063           % Some initialization of data
2064           if mdo.DCMODEL
2065             mpc.bus(:, VM) = 1;
2066           end
2067           % Injections and shadow prices
2068           mpc.gen(:, PG) = baseMVA * mdo.QP.x(vv.i1.Pg(t,j,k):vv.iN.Pg(t,j,k));
2069           %% need to update Qg for loads consistent w/constant power factor
2070           Pmin = mpc.gen(:, PMIN);
2071           Qmin = mpc.gen(:, QMIN);
2072           Qmax = mpc.gen(:, QMAX);
2073           ivl = find( isload(mpc.gen) & (Qmin ~= 0 | Qmax ~= 0) );
2074           Qlim = (Qmin(ivl) == 0) .* Qmax(ivl) + (Qmax(ivl) == 0) .* Qmin(ivl);
2075           mpc.gen(ivl, QG) = mpc.gen(ivl, PG) .* Qlim ./ Pmin(ivl);
2076           if mdo.DCMODEL
2077             %% bus angles
2078             mpc.bus(:, VA) = (180/pi) * mdo.QP.x(vv.i1.Va(t,j,k):vv.iN.Va(t,j,k));
2079           
2080             %% nodal prices
2081             price = (mdo.QP.lambda.mu_u(ll.i1.Pmis(t,j,k):ll.iN.Pmis(t,j,k))-mdo.QP.lambda.mu_l(ll.i1.Pmis(t,j,k):ll.iN.Pmis(t,j,k))) / baseMVA;
2082             mpc.bus(:, LAM_P) = price;
2083           
2084             %% line flows and line limit shadow prices
2085             mpc.branch(:, PF) = 0;
2086             mpc.branch(:, QF) = 0;
2087             mpc.branch(:, PT) = 0;
2088             mpc.branch(:, QT) = 0;
2089             mpc.branch(:, MU_SF) = 0;
2090             mpc.branch(:, MU_ST) = 0;
2091             ion = find(mpc.branch(:, BR_STATUS));
2092             rows = ll.i1.Pf(t,j,k):ll.iN.Pf(t,j,k);
2093             cols = vv.i1.Va(t,j,k):vv.iN.Va(t,j,k);
2094             lf = baseMVA * (mdo.QP.A(rows,cols) * mdo.QP.x(cols) + mdo.flow(t,j,k).PLsh);
2095             mpc.branch(ion, PF) = lf;
2096             mpc.branch(ion, PT) = -lf;
2097             mpc.branch(ion, MU_SF) = mdo.QP.lambda.mu_u(rows) / baseMVA;
2098             mpc.branch(ion, MU_ST) = mdo.QP.lambda.mu_l(rows) / baseMVA;
2099           else
2100             %% system price
2101             price = (mdo.QP.lambda.mu_l(ll.i1.Pmis(t,j,k):ll.iN.Pmis(t,j,k))-mdo.QP.lambda.mu_u(ll.i1.Pmis(t,j,k):ll.iN.Pmis(t,j,k))) / baseMVA;
2102             mpc.bus(:, LAM_P) = price;
2103           end
2104           if UC
2105             % igenon does not contain gens ousted because of a contingency or
2106             % a forced-off UC.CommitKey
2107             igenon = find(mpc.gen(:, GEN_STATUS));
2108             u = mdo.QP.x(vv.i1.u(t):vv.iN.u(t));
2109             mpc.gen(igenon, GEN_STATUS) = u(igenon);
2110             gs = mpc.gen(igenon, GEN_STATUS) > 0; % gen status
2111             mpc.gen(:, MU_PMAX) = 0;
2112             mpc.gen(:, MU_PMIN) = 0;
2113             mpc.gen(igenon, MU_PMAX) = gs .* ...
2114                     mdo.QP.lambda.mu_u(ll.i1.uPmax(t,j,k):ll.iN.uPmax(t,j,k)) / baseMVA;
2115             mpc.gen(igenon, MU_PMIN) = gs .* ...
2116                     mdo.QP.lambda.mu_u(ll.i1.uPmin(t,j,k):ll.iN.uPmin(t,j,k)) / baseMVA;
2117             if mdo.QCoordination
2118               mpc.gen(:, MU_QMAX) = 0;
2119               mpc.gen(:, MU_QMIN) = 0;
2120               mpc.gen(igenon, MU_QMAX) = gs .* ...
2121                       mdo.QP.lambda.mu_u(ll.i1.uQmax(t,j,k):ll.iN.uQmax(t,j,k)) / baseMVA;
2122               mpc.gen(igenon, MU_QMIN) = gs .* ...
2123                       mdo.QP.lambda.mu_u(ll.i1.uQmin(t,j,k):ll.iN.uQmin(t,j,k)) / baseMVA;
2124             end
2125           else
2126             gs = mpc.gen(:, GEN_STATUS) > 0;      % gen status
2127             mpc.gen(:, MU_PMAX) = gs .* ...
2128                     mdo.QP.lambda.upper(vv.i1.Pg(t,j,k):vv.iN.Pg(t,j,k)) / baseMVA;
2129             mpc.gen(:, MU_PMIN) = gs .* ...
2130                     mdo.QP.lambda.lower(vv.i1.Pg(t,j,k):vv.iN.Pg(t,j,k)) / baseMVA;
2131             if mdo.QCoordination
2132               mpc.gen(:, MU_QMAX) = gs .* ...
2133                       mdo.QP.lambda.upper(vv.i1.Qg(t,j,k):vv.iN.Qg(t,j,k)) / baseMVA;
2134               mpc.gen(:, MU_QMIN) = gs .* ...
2135                       mdo.QP.lambda.lower(vv.i1.Qg(t,j,k):vv.iN.Qg(t,j,k)) / baseMVA;
2136             end
2137           end
2138           if mdi.IncludeFixedReserves
2139             z = zeros(ng, 1);
2140             r = mdo.FixedReserves(t,j,k);
2141             r.R   = z;
2142             r.prc = z;
2143             r.mu = struct('l', z, 'u', z, 'Pmax', z);
2144             r.totalcost = sum(om.eval_quad_cost(mdo.QP.x, 'Rcost', {t,j,k}));
2145             r.R(r.igr) = mdo.QP.x(vv.i1.R(t,j,k):vv.iN.R(t,j,k)) * baseMVA;
2146             for gg = r.igr
2147               iz = find(r.zones(:, gg));
2148               kk = ll.i1.Rreq(t,j,k):ll.iN.Rreq(t,j,k);
2149               r.prc(gg) = sum(mdo.QP.lambda.mu_l(kk(iz))) / baseMVA;
2150             end
2151             r.mu.l(r.igr)    = mdo.QP.lambda.lower(vv.i1.R(t,j,k):vv.iN.R(t,j,k)) / baseMVA;
2152             r.mu.u(r.igr)    = mdo.QP.lambda.upper(vv.i1.R(t,j,k):vv.iN.R(t,j,k)) / baseMVA;
2153             r.mu.Pmax(r.igr) = mdo.QP.lambda.mu_u(ll.i1.Pg_plus_R(t,j,k):ll.iN.Pg_plus_R(t,j,k)) / baseMVA;
2154             mpc.reserves = r;
2155           end
2156           mdo.flow(t,j,k).mpc = mpc;     %% stash modified mpc in output struct
2157         end
2158       end
2159       % Contract, contingency reserves, energy limits
2160       mdo.results.Pc(:,t)  = baseMVA * mdo.QP.x(vv.i1.Pc(t):vv.iN.Pc(t));
2161       mdo.results.Rpp(:,t) = baseMVA * mdo.QP.x(vv.i1.Rpp(t):vv.iN.Rpp(t));
2162       mdo.results.Rpm(:,t) = baseMVA * mdo.QP.x(vv.i1.Rpm(t):vv.iN.Rpm(t));
2163       if ns
2164         mdo.results.Sm(:,t)  = baseMVA * mdo.QP.x(vv.i1.Sm(t):vv.iN.Sm(t));
2165         mdo.results.Sp(:,t)  = baseMVA * mdo.QP.x(vv.i1.Sp(t):vv.iN.Sp(t));
2166       end
2167     end
2168     % Ramping reserves
2169     for t = 1:mdo.idx.ntramp
2170       mdo.results.Rrp(:,t) = baseMVA * mdo.QP.x(vv.i1.Rrp(t):vv.iN.Rrp(t));
2171       mdo.results.Rrm(:,t) = baseMVA * mdo.QP.x(vv.i1.Rrm(t):vv.iN.Rrm(t));
2172     end
2173     % Expected energy prices for generators, per generator and per period,
2174     % both absolute and conditional on making it to that period
2175     mdo.results.GenPrices = zeros(ng, nt);
2176     mdo.results.CondGenPrices = zeros(ng, nt);
2177     for t = 1:nt
2178       pp = zeros(ng,1);
2179       for j = 1:mdo.idx.nj(t)
2180         for k = 1:mdo.idx.nc(t,j)+1
2181           pp = pp + mdo.flow(t,j,k).mpc.bus(mdo.flow(t,j,k).mpc.gen(:,GEN_BUS), LAM_P);
2182         end
2183       end
2184       mdo.results.GenPrices(:,t) = pp;
2185       mdo.results.CondGenPrices(:, t) = pp / mdo.StepProb(t);
2186     end
2187     % Obtain contingency reserve prices, per generator and period
2188     mdo.results.RppPrices = zeros(ng, nt);
2189     mdo.results.RpmPrices = zeros(ng, nt);
2190     for t = 1:nt
2191       mdo.results.RppPrices(:, t) = mdo.QP.lambda.lower(vv.i1.Rpp(t):vv.iN.Rpp(t)) / baseMVA;
2192       mdo.results.RpmPrices(:, t) = mdo.QP.lambda.lower(vv.i1.Rpm(t):vv.iN.Rpm(t)) / baseMVA;
2193       for j = 1:mdi.idx.nj(t);
2194         for k = 1:mdi.idx.nc(t,j)+1
2195           ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
2196           mdo.results.RppPrices(ii, t) = mdo.results.RppPrices(ii, t) + mdo.QP.lambda.mu_l(ll.i1.dPpRp(t,j,k):ll.iN.dPpRp(t,j,k)) / baseMVA;
2197           mdo.results.RpmPrices(ii, t) = mdo.results.RpmPrices(ii, t) + mdo.QP.lambda.mu_l(ll.i1.dPmRm(t,j,k):ll.iN.dPmRm(t,j,k)) / baseMVA;
2198         end
2199       end
2200     end
2201     % Obtain ramping reserve prices, per generator and period
2202     mdo.results.RrpPrices = zeros(ng, mdo.idx.ntramp);
2203     mdo.results.RrmPrices = zeros(ng, mdo.idx.ntramp);
2204     % First, 1:nt-1
2205     for t = 1:nt-1
2206       for j1 = 1:mdo.idx.nj(t)
2207         for j2 = 1:mdo.idx.nj(t+1)
2208           if mdi.tstep(t+1).TransMask(j2,j1)
2209             mdo.results.RrpPrices(:, t) = mdo.results.RrpPrices(:, t) + mdo.QP.lambda.mu_l(ll.i1.Rrp(t,j1,j2):ll.iN.Rrp(t,j1,j2)) / baseMVA;
2210             mdo.results.RrmPrices(:, t) = mdo.results.RrmPrices(:, t) + mdo.QP.lambda.mu_l(ll.i1.Rrm(t,j1,j2):ll.iN.Rrm(t,j1,j2)) / baseMVA;
2211           end
2212         end
2213       end
2214     end
2215     % then last period only if specified for with terminal state
2216     if ~mdo.OpenEnded
2217       for j1 = 1:mdo.idx.nj(nt)
2218         mdo.results.RrpPrices(:, nt) = mdo.results.RrpPrices(:, nt) + mdo.QP.lambda.mu_l(ll.i1.Rrp(nt,j1,1):ll.iN.Rrp(nt,j1,1)) / baseMVA;
2219         mdo.results.RrmPrices(:, nt) = mdo.results.RrmPrices(:, nt) + mdo.QP.lambda.mu_l(ll.i1.Rrm(nt,j1,1):ll.iN.Rrm(nt,j1,1)) / baseMVA;
2220       end
2221     end
2222     % Expected wear and tear costs per gen and period
2223     mdo.results.ExpectedRampCost = zeros(ng, mdo.idx.ntramp+1);
2224     % First do first period wrt to InitialPg.
2225     for j = 1:mdi.idx.nj(1)
2226       w = mdo.tstep(1).TransMat(j,1); % the probability of going from initial state to jth
2227       mdo.results.ExpectedRampCost(:, 1) = mdo.results.ExpectedRampCost(:, 1) ...
2228           + 0.5 * w * mdo.RampWearCostCoeff(:,1) .* (mdo.flow(1,j,1).mpc.gen(:,PG) - mdo.InitialPg).^2;
2229     end
2230     % Then the remaining periods
2231     for t = 2:nt
2232       for j2 = 1:mdo.idx.nj(t)
2233         for j1 = 1:mdo.idx.nj(t-1)
2234           w = mdo.tstep(t).TransMat(j2,j1) * mdo.CostWeights(1, j1, t-1);
2235           mdo.results.ExpectedRampCost(:, t) = mdo.results.ExpectedRampCost(:, t) ...
2236               + 0.5 * w * mdo.RampWearCostCoeff(:,t) .* (mdo.flow(t,j2,1).mpc.gen(:,PG) - mdo.flow(t-1,j1,1).mpc.gen(:,PG)) .^2;
2237         end
2238       end
2239     end
2240     % Finally, if there is a terminal state problem, apply cost to
2241     if ~mdo.OpenEnded
2242       for j = 1:mdi.idx.nj(nt)
2243         w = mdi.tstep(t+1).TransMat(1, j) * mdi.CostWeights(1, j, nt);
2244         mdo.results.ExpectedRampCost(:, nt+1) = 0.5 * w * mdo.RampWearCostCoeff(:,nt+1) .* (mdo.TerminalPg - mdo.flow(nt,j,1).mpc.gen(:,PG)) .^2;
2245       end
2246     end
2247     % Compute expected dispatch, conditional on making it to the
2248     % corresponding period
2249     mdo.results.ExpectedDispatch = zeros(ng, nt);
2250     for t = 1:nt
2251       pp = sum(mdo.CostWeights(1,1:mdo.idx.nj(t),t)');    % gamma(t+1)
2252       for j = 1:mdo.idx.nj(t)
2253         mdo.results.ExpectedDispatch(:,t) = mdo.results.ExpectedDispatch(:,t) + ...
2254               mdo.CostWeights(1,j,t)/pp * mdo.flow(t,j,1).mpc.gen(:,PG);
2255       end
2256     end
2257     % If Cyclic storage, pull InitialStorage value out of x
2258     if ns && mdo.Storage.ForceCyclicStorage
2259       mdo.Storage.InitialStorage = baseMVA * mdo.QP.x(vv.i1.S0:vv.iN.S0);
2260     end
2261     % Compute expected storage state trajectory
2262     mdo.Storage.ExpectedStorageState = zeros(ns,nt);
2263     if ns
2264       for t = 1:nt
2265         pp = sum(mdo.CostWeights(1,1:mdo.idx.nj(t),t)');    %% gamma(t+1)
2266         for j = 1:mdo.idx.nj(t)
2267           Lfj = mdo.tstep(t).Lf( j:mdo.idx.nj(t):(ns-1)*mdo.idx.nj(t)+j, :);
2268           Ngj = mdo.tstep(t).Ng( j:mdo.idx.nj(t):(ns-1)*mdo.idx.nj(t)+j, :);
2269           Nhj = mdo.tstep(t).Nh( j:mdo.idx.nj(t):(ns-1)*mdo.idx.nj(t)+j, :);
2270           mdo.Storage.ExpectedStorageState(:,t) = ...
2271               mdo.Storage.ExpectedStorageState(:,t) + ...
2272                   baseMVA * mdo.CostWeights(1,j,t)/pp * ...
2273                   ( Lfj * mdo.Storage.InitialStorage/baseMVA + ...
2274                     (Ngj + Nhj) * mdo.QP.x );
2275         end
2276       end
2277       mdo.Storage.ExpectedStorageDispatch = ...
2278           mdo.results.ExpectedDispatch(mdo.Storage.UnitIdx, :);
2279     end
2280     % If there is a dynamical system, extract the state vectors and outputs
2281     % from the solution
2282     if ntds
2283       if nzds
2284         mdo.results.Z = zeros(nzds, ntds);
2285         for t = 1:ntds
2286           mdo.results.Z(:,t) = mdo.QP.x(vv.i1.Z(t):vv.iN.Z(t));
2287         end
2288       end
2289       mdo.results.Y = zeros(nyds, ntds);
2290       if nyds
2291         for t = 1:ntds
2292           mdo.results.Y(:, t) = ...
2293                   mdo.QP.A(ll.i1.DSy(t):ll.iN.DSy(t), :) * mdo.QP.x;
2294         end
2295       end
2296     end
2297     mdo.results.f = mdo.QP.f;
2298   end   % if success
2299 end     % if mpopt.most.solve_model
2300 
2301 mdo.results.success = success;
2302 mdo.results.SetupTime = et_setup;
2303 mdo.results.SolveTime = toc(t0);
2304 
2305 if verbose
2306   fprintf('- MOST: Done.\n\n');
2307 end