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 OPT_MODEL/SOLVE
                   (i.e. MIQPS_MASTER or QPS_MASTER) for details
           most.skip_prices (0) - skip price computation stage for mixed
                   integer problems, see MIQPS_MASTER 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 MIQPS_MASTER 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 OPT_MODEL/SOLVE
0048 %                   (i.e. MIQPS_MASTER or QPS_MASTER) for details
0049 %           most.skip_prices (0) - skip price computation stage for mixed
0050 %                   integer problems, see MIQPS_MASTER 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 MIQPS_MASTER 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-2020, 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-2020 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   if ns
0712     for t = 1:nt
0713       for j = 1:mdi.idx.nj(t)
0714         for k = 1:mdi.idx.nc(t,j)+1
0715           om.add_var('Psc', {t,j,k}, ns, [], [], zeros(ns,1));
0716         end
0717       end
0718     end
0719     for t = 1:nt
0720       for j = 1:mdi.idx.nj(t)
0721         for k = 1:mdi.idx.nc(t,j)+1
0722           om.add_var('Psd', {t,j,k}, ns, [], zeros(ns,1), []);
0723         end
0724       end
0725     end
0726   end
0727   % Continue with storage upper and lower bounds, one for each time period
0728   % and unit
0729   om.init_indexed_name('var', 'Sp', {nt});
0730   om.init_indexed_name('var', 'Sm', {nt});
0731   if ns
0732     for t = 1:nt
0733       om.add_var('Sp', {t}, ns, [], [], MaxStorageLevel(:,t)/baseMVA);
0734     end
0735     for t = 1:nt
0736       om.add_var('Sm', {t}, ns, [], MinStorageLevel(:,t)/baseMVA, []);
0737     end
0738   end
0739   % Possible initial storage quantities when using cyclic storage dispatch
0740   % so that initial storage = expected terminal storage is a constraint
0741   if ns && mdi.Storage.ForceCyclicStorage
0742     om.add_var('S0', ns, [], ...
0743         mdi.Storage.InitialStorageLowerBound / baseMVA, ...
0744         mdi.Storage.InitialStorageUpperBound / baseMVA);
0745   end
0746   % If there is a dynamical system with non-null state vector,
0747   % add those states here
0748   if nzds
0749     om.init_indexed_name('var', 'Z', {ntds});
0750     for t = 1:ntds
0751       if t == 1
0752         zmin = mdi.z1;
0753         zmax = mdi.z1;
0754       else
0755         zmin = mdi.dstep(t).zmin;
0756         zmax = mdi.dstep(t).zmax;
0757       end
0758       z0 = (zmax - zmin) / 2;
0759       om.add_var('Z', {t}, nzds, z0, zmin, zmax);
0760     end
0761   end
0762   % Now the integer variables; u variables mean on/off status
0763   if UC
0764     om.init_indexed_name('var', 'u', {nt});
0765     om.init_indexed_name('var', 'v', {nt});
0766     om.init_indexed_name('var', 'w', {nt});
0767     vt0 = char('B' * ones(1, ng));  % default variable type for u is binary
0768     for t = 1:nt
0769       umin = zeros(ng, 1);
0770       umax = ones(ng, 1);
0771       % min up/down restrictions on u
0772       if ~mdi.UC.CyclicCommitment
0773         % min up time has not passed yet since startup occured, force ON
0774         umin( (mdi.UC.InitialState > 0) & ...
0775               (t+mdi.UC.InitialState-mdi.UC.MinUp <= 0) ) = 1;
0776         % min down time has not passed yet since shutdown occured, force OFF
0777         umax( (mdi.UC.InitialState < 0) & ...
0778               (t-mdi.UC.InitialState-mdi.UC.MinDown <= 0) ) = 0;
0779       end
0780       % set limits for units forced ON or forced OFF
0781       iON = find(mdi.UC.CommitKey(:,t) == 2);
0782       iOFF = find(mdi.UC.CommitKey(:,t) == -1);
0783       umin(iON)  = 1;
0784       umax(iOFF) = 0;
0785 
0786       % set variable types
0787       vt = vt0;                 % initialize all variable types to binary
0788       vt(umin == umax) = 'C';   % make continuous for those that are fixed
0789       
0790       om.add_var('u', {t}, ng, zeros(ng, 1), umin, umax, vt);
0791     end
0792     % v variables mean startup events
0793     for t = 1:nt
0794       om.add_var('v', {t}, ng, zeros(ng, 1), zeros(ng, 1), ones(ng, 1));
0795     end
0796     % w variables mean shutdown events
0797     for t = 1:nt
0798       om.add_var('w', {t}, ng, zeros(ng, 1), zeros(ng, 1), ones(ng, 1));
0799     end
0800   end
0801   % An external program may be using coordination with AC flows, and in
0802   % that case we need corresponding Qg variables, whose only function is to
0803   % be constrained to zero if the commitment decision asks for that.
0804   if mdi.QCoordination
0805     om.init_indexed_name('var', 'Qg', {nt, nj_max, nc_max+1});
0806     for t = 1:nt
0807       for j = 1:mdi.idx.nj(t)
0808         for k = 1:mdi.idx.nc(t,j)+1
0809           mpc = mdi.flow(t,j,k).mpc;
0810           genmask = mpc.gen(:,GEN_STATUS) > 0;
0811           q0 = genmask .* mpc.gen(:, QG) / baseMVA;
0812           if UC     % relax bounds here, enforced by uQmax, uQmin constraints
0813             qmin = genmask .* (min(mpc.gen(:, QMIN) / baseMVA, 0) - 1);
0814             qmax = genmask .* (max(mpc.gen(:, QMAX) / baseMVA, 0) + 1);
0815           else      % enforce bounds here, subject to flow's GEN_STATUS
0816             qmin = genmask .* mpc.gen(:, QMIN) / baseMVA;
0817             qmax = genmask .* mpc.gen(:, QMAX) / baseMVA;
0818           end
0819           om.add_var('Qg', {t,j,k}, ng, q0, qmin, qmax);
0820         end
0821       end
0822     end
0823   end
0824 
0825   %% handing of user-defined variables would go here
0826 
0827   nvars = om.getN('var');
0828   mdi.idx.nvars = nvars;
0829 
0830   % Construct mechanism to keep track of expected storage states. This is
0831   % used for building some constraints in some types of problems, most
0832   % notably when there is a constraint on the terminal expected storage,
0833   % but it is also needed when there is a value associated with leftover
0834   % storage. Unfortunately this requires the construction of a substantial
0835   % mechanism for computing the expected terminal storage at any end-point
0836   % of the transition tree, be it a contingency or the terminal of the central
0837   % path at the end of the horizon.  Let SF(t) be the terminal storage value
0838   % at the end of the t-th period, assuming that we make it there without
0839   % contingencies, and let SI(t) be the initial amount of storage at that
0840   % same period. Then
0841   %     SI(t) = D(t) * SF(t-1).
0842   %     SF(t) = B1(t) * SI(t) + B2*[G(t)*x + H(t)]*x, and
0843   % Here, D(t) is created from the probability transition matrix, restricted
0844   % to transitions from base cases at t-1 to base cases at t. This allows
0845   % us to write a recursion. If the general form for SI(t),SF(t) is
0846   %     SI(t) = Li(t)*S0 + Mg(t)*x + Mh(t)*x,
0847   %     SF(t) = Lf(t)*S0 + Ng(t)*x + Nh(t)*x,
0848   % it turns out that the recursion is
0849   %     L(1) = D(1); Mg(1) = Mh(1) = 0; Ng(1) = G(1); Nh(1) = H(1);
0850   %     for t=2:nt
0851   %         Li(t) = D(t)*Lf(t-1) = D(t)*B1(t-1)*Li(t-1);
0852   %         Lf(t) = B1(t)*Li(t) = B1(t)*D(t)*Lf(t-1);
0853   %         Mg(t) = D(t)*Ng(t-1);
0854   %         Mh(t) = D(t)*Nh(t-1);
0855   %         Ng(t) = B1(t)*Mg(t) + B2(t)*G(t);
0856   %         Nh(t) = B1(t)*Mh(t) + B2(t)*H(t);
0857   %     end
0858   %
0859   % If SI,SF are organized first by blocks for each storage unit and within
0860   % the blocks by scenario, then the D matrix is simply made up by
0861   % repeating the str.tstep(t).TransMat matrix ns times in the diagonal
0862   % and then the columns weighted by the probabilities of the basecases at
0863   % t-1 given that we remained in basecases, and the rows are weighted by
0864   % the inverse of the probabilites of the scenarios at the beginning of t.
0865   % D(1) is special though, where each block is an nj(1) x 1 vector of ones.
0866   % The B matrices are formed by stacking appropriately sized diagonal
0867   % matrices for each storage unit along the diagonal, where each component
0868   % is simply the i-th element of beta times an identity matrix.
0869   if verbose
0870     fprintf('- Building expected storage-tracking mechanism.\n');
0871   end
0872   if ns
0873     % The following code assumes that no more variables will be added
0874     vv = om.get_idx();
0875     for t = 1:nt
0876       nsxnjt = ns*mdi.idx.nj(t);
0877       % Form G(t), H(t), B1(t), B2(t)
0878       G = sparse(nsxnjt, nvars);
0879       H = sparse(nsxnjt, nvars);
0880       B1 = sparse(nsxnjt, nsxnjt);
0881       B2 = sparse(nsxnjt, nsxnjt);
0882       for j = 1:mdi.idx.nj(t)
0883         ii  = ((1:ns)'-1)*mdi.idx.nj(t)+j;
0884         jj1 = (vv.i1.Psc(t,j,1):vv.iN.Psc(t,j,1))';
0885         jj2 = (vv.i1.Psd(t,j,1):vv.iN.Psd(t,j,1))';
0886         G = G + sparse(ii, jj1, -mdi.Delta_T  *  InEff(:,t), nsxnjt, nvars);
0887         H = H + sparse(ii, jj2, -mdi.Delta_T ./ OutEff(:,t), nsxnjt, nvars);
0888         B1 = B1 + sparse(ii, ii, beta1(:,t), nsxnjt, nsxnjt);
0889         B2 = B2 + sparse(ii, ii, beta2(:,t), nsxnjt, nsxnjt);
0890       end
0891       if t == 1
0892         % form Li, Lf, Mg, Mh, Ng, Nh, B1, B2 for t == 1
0893         jlist = [];
0894         for i=1:ns
0895           jlist = [ jlist; i*ones(mdi.idx.nj(t),1) ];
0896         end
0897         mdi.tstep(t).Li  = sparse((1:nsxnjt)', jlist, 1, nsxnjt, ns);
0898         mdi.tstep(t).Lf  = B1 * mdi.tstep(t).Li;
0899         mdi.tstep(t).Mg  = sparse(nsxnjt, nvars);   % Initial one is all zeros
0900         mdi.tstep(t).Mh  = sparse(nsxnjt, nvars);   % Initial one is all zeros
0901         mdi.tstep(t).Ng  = B2 * G;
0902         mdi.tstep(t).Nh  = B2 * H;
0903       else
0904         % Form D(t)
0905         D = sparse(nsxnjt, ns*mdi.idx.nj(t-1));
0906         p1 = mdi.CostWeights(1,1:mdi.idx.nj(t-1),t-1)';
0907         p1 = p1 / sum(p1);      % sigma(t)
0908         p2 = mdi.tstep(t).TransMat * p1;
0909         Di = spdiags(1./p2, 0, mdi.idx.nj(t), mdi.idx.nj(t)) * ...
0910                 sparse(mdi.tstep(t).TransMat) * ...
0911                 spdiags(p1, 0, mdi.idx.nj(t-1), mdi.idx.nj(t-1));
0912         for i = 1:ns
0913           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;
0914         end
0915         % Apply recursion, form Li, Lf, Mg, Mh, Ng, Nh
0916         mdi.tstep(t).Li = D  * mdi.tstep(t-1).Lf;
0917         mdi.tstep(t).Lf = B1 * mdi.tstep(t).Li;
0918         mdi.tstep(t).Mg = D * mdi.tstep(t-1).Ng;
0919         mdi.tstep(t).Mh = D * mdi.tstep(t-1).Nh;
0920         mdi.tstep(t).Ng = B1 * mdi.tstep(t).Mg + B2 * G;
0921         mdi.tstep(t).Nh = B1 * mdi.tstep(t).Mh + B2 * H;
0922       end
0923       mdi.tstep(t).G = G;
0924       mdi.tstep(t).H = H;
0925     end
0926   end
0927 
0928   % Now for the constraint indexing and creation.
0929   if verbose
0930     fprintf('- Building constraint submatrices.\n');
0931   end
0932   baseMVA = mdi.mpc.baseMVA;
0933   om.init_indexed_name('lin', 'Pmis', {nt, nj_max, nc_max+1});
0934   if mdi.DCMODEL
0935     % Construct all load flow equations using a DC flow model
0936     if verbose
0937       fprintf('  - Building DC flow constraints.\n');
0938     end
0939     om.init_indexed_name('lin', 'Pf', {nt, nj_max, nc_max+1});
0940     for t = 1:nt
0941       for j = 1:mdi.idx.nj(t)
0942         for k = 1:mdi.idx.nc(t,j)+1
0943           % First the flow constraints
0944           mpc = mdi.flow(t,j,k).mpc;
0945           ion = find(mpc.branch(:, BR_STATUS));
0946           [Bdc, Bl, Psh, PLsh] = makeBdc(baseMVA, mpc.bus, mpc.branch(ion,:));
0947           mdi.flow(t,j,k).PLsh = PLsh;     %% save for computing flows later
0948           negCg = sparse(mpc.gen(:,GEN_BUS), (1:ng)', -1, ...
0949                         mdi.idx.nb(t,j,k), ng);
0950           A = [Bdc negCg];
0951           b = -(mpc.bus(:,PD)+mpc.bus(:,GS))/baseMVA-Psh;
0952           vs = struct('name', {'Va', 'Pg'}, 'idx', {{t,j,k}, {t,j,k}});
0953           om.add_lin_constraint('Pmis', {t,j,k}, A, b, b, vs);
0954           % Then the thermal limits
0955           tmp = mpc.branch(ion,RATE_A)/baseMVA;
0956           iuncon = find(~tmp);
0957           tmp(iuncon) = Inf(size(iuncon));
0958           vs = struct('name', {'Va'}, 'idx', {{t,j,k}});
0959           om.add_lin_constraint('Pf', {t,j,k}, Bl, -tmp-PLsh, tmp-PLsh, vs);
0960         end
0961       end
0962     end
0963   else
0964     if verbose
0965       fprintf('  - Building load balance constraints.\n');
0966     end
0967     % Set simple generation - demand = 0 equations, one for each flow
0968     for t = 1:nt
0969       for j = 1:mdi.idx.nj(t)
0970         for k = 1:mdi.idx.nc(t,j)+1
0971           mpc = mdi.flow(t,j,k).mpc;
0972           A = sparse(ones(1, ng));
0973           b = 1.0*sum(mpc.bus(:, PD)+mpc.bus(:,GS))/baseMVA;
0974           vs = struct('name', {'Pg'}, 'idx', {{t,j,k}});
0975           om.add_lin_constraint('Pmis', {t,j,k}, A, b, b, vs);
0976         end
0977       end
0978     end
0979   end
0980   if mdi.IncludeFixedReserves
0981     if verbose
0982       fprintf('  - Building fixed zonal reserve constraints.\n');
0983     end
0984     om.init_indexed_name('lin', 'Pg_plus_R', {nt, nj_max, nc_max+1});
0985     om.init_indexed_name('lin', 'Rreq', {nt, nj_max, nc_max+1});
0986     for t = 1:nt
0987       for j = 1:mdi.idx.nj(t)
0988         for k = 1:mdi.idx.nc(t,j)+1
0989           % First the flow constraints
0990           mpc = mdi.flow(t,j,k).mpc;
0991           r = mdi.FixedReserves(t,j,k);
0992           ngr = length(r.igr);
0993           I = speye(ngr);
0994           Ar = sparse(1:ngr, r.igr, 1, ngr, ng);
0995           if UC
0996             A = [Ar I  ...
0997                   sparse(1:ngr, r.igr, -mpc.gen(r.igr, PMAX) / baseMVA, ngr, ng)];
0998             u = zeros(ngr, 1);
0999             vs = struct('name', {'Pg', 'R', 'u'}, 'idx', {{t,j,k}, {t,j,k}, {t}});
1000           else
1001             A = [Ar I];
1002             u = mpc.gen(r.igr, PMAX) / baseMVA;
1003             vs = struct('name', {'Pg', 'R'}, 'idx', {{t,j,k}, {t,j,k}});
1004           end
1005           om.add_lin_constraint('Pg_plus_R', {t,j,k}, A, [], u, vs);
1006           A = r.zones(:, r.igr);
1007           l = r.req / mpc.baseMVA;
1008           vs = struct('name', {'R'}, 'idx', {{t,j,k}});
1009           om.add_lin_constraint('Rreq', {t,j,k}, A, l, [], vs);
1010         end
1011       end
1012     end
1013   end
1014   
1015   % Set relationships between generator injections and charge/discharge
1016   % variables (-pg + psc + psd = 0)
1017   if verbose && ~isempty(mdi.Storage.UnitIdx)
1018     fprintf('  - Splitting storage injections into charge/discharge.\n');
1019   end
1020   om.init_indexed_name('lin', 'Ps', {nt, nj_max, nc_max+1});
1021   for t = 1:nt
1022     for j = 1:mdi.idx.nj(t)
1023       for k = 1:mdi.idx.nc(t,j)+1
1024         A = [sparse((1:ns)', mdi.Storage.UnitIdx, -1, ns, ng) Ins Ins];
1025         b = zeros(ns, 1);
1026         vs = struct('name', {'Pg', 'Psc', 'Psd'}, 'idx', {{t,j,k}, {t,j,k}, {t,j,k}});
1027         om.add_lin_constraint('Ps', {t,j,k}, A, b, b, vs);
1028       end
1029     end
1030   end
1031   
1032   % Construct y-variable restrictions on piecewise-linear costs. Note that
1033   % the restriction lines are computed using the full non-scaled cost in
1034   % gencost; any weighting of the cost must be then specified later in the
1035   % cost coefficients hitting the y variables (not 1 anymore).  Do it for
1036   % a complete c3sopf cell taking advantage of the fact that all p
1037   % injections are contiguous, as are all y variables for a c3sopf cell.
1038   % Also, note that makeAy assumes that every gen is online, which is
1039   % consistent with our formulation for unit commitment
1040   if verbose
1041     fprintf('  - Building CCV constraints for piecewise-linear costs.\n');
1042   end
1043   om.init_indexed_name('lin', 'ycon', {nt, nj_max, nc_max+1});
1044   for t = 1:nt
1045     for j = 1:mdi.idx.nj(t)
1046       for k = 1:mdi.idx.nc(t,j)+1
1047         mpc = mdi.flow(t,j,k).mpc;
1048         [A, u] = makeAy(baseMVA, ng, mpc.gencost, 1, [], ng+1);
1049         vs = struct('name', {'Pg', 'y'}, 'idx', {{t,j,k}, {t,j,k}});
1050         om.add_lin_constraint('ycon', {t,j,k}, A, [], u, vs);
1051       end
1052     end
1053   end
1054 
1055   % The actual deviations from base flow must not exceed physical ramp rates
1056   % we'll get negative multiplier for right bound, fix when picking up
1057   % lambdas.
1058   % At issue: generators ousted in a contingency clearly violate this
1059   % transition; do not include constraints for these or for generators
1060   % whose commitment key is -1; either of these two possibilities will
1061   % result in a GEN_STATUS of 0, so we use that as the indicator.
1062   if verbose
1063     fprintf('  - Building contingency reserve constraints.\n');
1064   end
1065   om.init_indexed_name('lin', 'rampcont', {nt, nj_max, nc_max+1});
1066   for t =1:nt
1067     for j = 1:mdi.idx.nj(t)
1068       for k = 2:mdi.idx.nc(t,j)+1
1069         ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
1070         ngtmp = length(ii);
1071         A = sparse((1:ngtmp)', ii, 1, ngtmp, ng);
1072         u = mdi.flow(t,j,k).mpc.gen(ii,RAMP_10)/baseMVA;
1073         vs = struct('name', {'Pg', 'Pg'}, 'idx', {{t,j,1}, {t,j,k}});
1074         om.add_lin_constraint('rampcont', {t,j,k}, [-A A], -u, u, vs);
1075       end
1076     end
1077   end
1078   % The usual alpha-controlled equality of P0 and Pc does not make
1079   % sense when having many scenarios and hence many P0's .  Ditch.
1080   %
1081   % Positive reserve variables are larger than all increment variables in
1082   % all scenarios and flows of a given time slice 0 <= rpp - dpp; these
1083   % are the ones that set the price of reserves. Include all units that are
1084   % potentially committed.
1085   om.init_indexed_name('lin', 'dPpRp', {nt, nj_max, nc_max+1});
1086   for t = 1:nt
1087     for j = 1:mdi.idx.nj(t);
1088       for k = 1:mdi.idx.nc(t,j)+1
1089         ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
1090         ngtmp = length(ii);
1091         A = sparse((1:ngtmp)', ii, 1, ngtmp, ng);
1092         l = zeros(ngtmp, 1);
1093         vs = struct('name', {'dPp', 'Rpp'}, 'idx', {{t,j,k}, {t}});
1094         om.add_lin_constraint('dPpRp', {t,j,k}, [-A A], l, [], vs);
1095       end
1096     end
1097   end
1098   % Negative reserve variables are larger than all decrement variables in
1099   % all scenarios and flows of a given time slice  0 <= rpm - dpm; these
1100   % are the ones that set the price of reserves. Include all units that are
1101   % potentially committed.
1102   om.init_indexed_name('lin', 'dPmRm', {nt, nj_max, nc_max+1});
1103   for t = 1:nt
1104     for j = 1:mdi.idx.nj(t);
1105       for k = 1:mdi.idx.nc(t,j)+1
1106         ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
1107         ngtmp = length(ii);
1108         A = sparse((1:ngtmp)', ii, 1, ngtmp, ng);
1109         l = zeros(ngtmp, 1);
1110         vs = struct('name', {'dPm', 'Rpm'}, 'idx', {{t,j,k}, {t}});
1111         om.add_lin_constraint('dPmRm', {t,j,k}, [-A A], l, [], vs);
1112       end
1113     end
1114   end
1115   % The difference between the injection and the contract
1116   % is equal to the inc minus the dec: Ptjk - Ptc = dPp - dPm
1117   % Include all units that are potentially committed.
1118   om.init_indexed_name('lin', 'dPdef', {nt, nj_max, nc_max+1});
1119   for t = 1:nt
1120     for j = 1:mdi.idx.nj(t)
1121       for k = 1:mdi.idx.nc(t,j)+1
1122         ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
1123         ngtmp = length(ii);
1124         A = sparse((1:ngtmp)', ii, 1, ngtmp, ng);
1125         b = zeros(ngtmp, 1);
1126         vs = struct('name', {'Pg', 'Pc', 'dPp', 'dPm'}, ...
1127                     'idx', {{t,j,k}, {t}, {t,j,k}, {t,j,k}});
1128         om.add_lin_constraint('dPdef', {t,j,k}, [A -A -A A], b, b, vs);
1129       end
1130     end
1131   end
1132   
1133   % Go on to load following ramping restrictions.  Note that these
1134   % restrictions apply even if there is a change in the commitment status
1135   % of the generator.
1136   %
1137   % First, bound upward ramping reserves from below by all base-case
1138   % ramping possibilities, 0 <= rrp(t) -p(t+1)(j2)0 + p(t)(j1)0.  A ramping
1139   % reserve is needed at time t to be able to change to the needed dispatch at
1140   % time t+1.  An initial ramping reserve (t=0) would be data, not a
1141   % variable, and while ramping transitions from t=0 to t=1 should be
1142   % enforced, we do not allocate the reserve for t = 0.
1143   % Note: in the event that some future reserves are already locked in, we may
1144   %       not want to start at t = 1
1145   if verbose
1146     fprintf('  - Building ramping transitions and reserve constraints.\n');
1147   end
1148   om.init_indexed_name('lin', 'Rrp', {nt, nj_max, nj_max});
1149   % First, do from t=1:nt-1, since the last one is different and may not
1150   % even exist depending on the type of horizon
1151   for t = 1:nt-1
1152     for j1 = 1:mdi.idx.nj(t)      % j1 is at time t
1153       for j2 = 1:mdi.idx.nj(t+1)  % j2 is at time t+1
1154         if mdi.tstep(t+1).TransMask(j2,j1)
1155           A = [Ing -Ing Ing];
1156           l = zeros(ng, 1);
1157           vs = struct('name', {'Pg', 'Pg', 'Rrp'}, ...
1158                       'idx', {{t,j1,1}, {t+1,j2,1}, {t}});
1159           om.add_lin_constraint('Rrp', {t,j1,j2}, A, l, [], vs);
1160         end
1161       end
1162     end
1163   end
1164   % Now, pay special attention to a possible last type of ramping
1165   % constraint. If the horizon involves a terminal value at t=nt+1, then
1166   % this must also be enforced; in this case, additional ramping
1167   % reserves procured for t=nt must be defined.  If this
1168   % condition does not apply, then these reserves are not needed.
1169   if ~mdi.OpenEnded
1170     % pterminal <= rrp(nt) + p(nt,j1,0)
1171     for j1 = 1:mdi.idx.nj(nt)
1172       A = [Ing Ing];
1173       l = mdi.TerminalPg/baseMVA;
1174       vs = struct('name', {'Pg', 'Rrp'}, ...
1175                   'idx', {{nt,j1,1}, {nt}});
1176       om.add_lin_constraint('Rrp', {nt,j1,1}, A, l, [], vs);
1177     end
1178   end
1179   % Now on to downward ramping reserves.
1180   % First, bound downward ramping reserves from below by all base-case
1181   % ramping possibilities, 0 <= rrm(t) + p(t+1)j20 - p(t)j10
1182   om.init_indexed_name('lin', 'Rrm', {nt, nj_max, nj_max});
1183   % First, do from t=1:nt-1, since the last one is different and may not
1184   % even exist depending on the type of horizon
1185   for t = 1:nt-1
1186     for j1 = 1:mdi.idx.nj(t)      % j1 is at time t
1187       for j2 = 1:mdi.idx.nj(t+1)  % j2 is at time t+1
1188         if mdi.tstep(t+1).TransMask(j2,j1)
1189           A = [-Ing Ing Ing];
1190           l = zeros(ng, 1);
1191           vs = struct('name', {'Pg', 'Pg', 'Rrm'}, ...
1192                       'idx', {{t,j1,1}, {t+1,j2,1}, {t}});
1193           om.add_lin_constraint('Rrm', {t,j1,j2}, A, l, [], vs);
1194         end
1195       end
1196     end
1197   end
1198   % Now, pay special attention to a possible last type of ramping
1199   % constraint. If the horizon involves a terminal value at t=nt+1, then
1200   % this must also be enforced; in this case, additional ramping
1201   % reserves procured for t=nt must be defined.  If this
1202   % condition does not apply, then these reserves are not needed.
1203   if ~mdi.OpenEnded
1204     % -pterminal <= rrm(nt) - p(nt,j1,0)
1205     for j1 = 1:mdi.idx.nj(nt)
1206       A = [-Ing Ing];
1207       l = -mdi.TerminalPg/baseMVA;
1208       vs = struct('name', {'Pg', 'Rrm'}, ...
1209                   'idx', {{nt,j1,1}, {nt}});
1210       om.add_lin_constraint('Rrm', {nt,j1,1}, A, l, [], vs);
1211     end
1212   end
1213   %
1214   %
1215   % Now for the storage restrictions.
1216   if ns
1217     if verbose
1218       fprintf('  - Building storage constraints.\n');
1219     end
1220     % First bound sm(t) based on sm(t-1), with sm(1) being bound by the initial
1221     % data; this is for base case trajectories only
1222     om.init_indexed_name('lin', 'Sm', {nt, nj_max});
1223     if mdi.Storage.ForceCyclicStorage
1224       % sm(1) - beta1*s0 + beta2*Delta_T*[eta_c*psc(1,j,0) + (1/eta_d)*psd(1,j,0)] <= 0
1225       for j = 1:mdi.idx.nj(1)
1226         A = [ diagBeta2EtaIn1 diagBeta2overEtaOut1 Ins -spdiags(beta1(:,1), 0, ns, ns)];
1227         u = zeros(ns, 1);
1228         vs = struct('name', {'Psc', 'Psd', 'Sm', 'S0'}, 'idx', {{1,j,1}, {1,j,1}, {1}, {}});
1229         om.add_lin_constraint('Sm', {1,j}, A, [], u, vs);
1230       end
1231     else
1232       % sm(1) + beta2*Delta_T*[eta_c*psc(1,j,0) + (1/eta_d)*psd(1,j,0)] <= beta1*(rho*InitialLB+(1-rho)*Initial)/baseMVA
1233       for j = 1:mdi.idx.nj(1)
1234         A = [ diagBeta2EtaIn1 diagBeta2overEtaOut1 Ins ];
1235         u = beta1(:,1) .* ( rho(:,t).*mdi.Storage.InitialStorageLowerBound + ...
1236                             (1-rho(:,t)).*mdi.Storage.InitialStorage ) / baseMVA;
1237         vs = struct('name', {'Psc', 'Psd', 'Sm'}, 'idx', {{1,j,1}, {1,j,1}, {1}});
1238         om.add_lin_constraint('Sm', {1,j}, A, [], u, vs);
1239       end
1240     end
1241     % Then the rest of the periods
1242     % 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
1243     % where s_I(t,j) = L_I(t,j) * s0 + (Mg(t,j)+Mh(t,j)) * x
1244     for t = 2:nt
1245       for j = 1:mdi.idx.nj(t)
1246         Mj = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :) + ...
1247              mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1248         Lij = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1249         diag1minusRhoBeta1 = spdiags((1-rho(:,t)) .* beta1(:,t), 0, ns, ns);
1250         A = sparse([1:ns,1:ns,1:ns,1:ns]', ...
1251                    [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)]', ...
1252                    [beta2EtaIn(:,t); beta2overEtaOut(:,t); -beta1(:,t).*rho(:,t); ones(ns,1)], ...
1253                    ns, nvars) ...
1254                 - diag1minusRhoBeta1 * Mj;
1255         if mdi.Storage.ForceCyclicStorage
1256           As0 = sparse(ns, nvars);
1257           As0(:, vv.i1.S0:vv.iN.S0) = -diag1minusRhoBeta1 * Lij;
1258           A = A + As0;
1259           u = zeros(ns, 1);
1260         else
1261           u = full(diag1minusRhoBeta1 * Lij * mdi.Storage.InitialStorage/baseMVA);
1262         end
1263         om.add_lin_constraint('Sm', {t,j}, A, [], u);
1264       end
1265     end
1266     % Do the same we did for sm(t) for sp(t). First the initial step ...
1267     om.init_indexed_name('lin', 'Sp', {nt, nj_max});
1268     if mdi.Storage.ForceCyclicStorage
1269       % -sp(1) + beta1*s0 - beta2*Delta_T*[eta_c*psc(1,j,0) + (1/eta_d)*psd(1,j,0)] <= 0
1270       for j = 1:mdi.idx.nj(1)
1271         A = [ -diagBeta2EtaIn1 -diagBeta2overEtaOut1 -Ins spdiags(beta1(:,1), 0, ns, ns) ];
1272         u = zeros(ns, 1);
1273         vs = struct('name', {'Psc', 'Psd', 'Sp', 'S0'}, 'idx', {{1,j,1}, {1,j,1}, {1}, {}});
1274         om.add_lin_constraint('Sp', {1,j}, A, [], u, vs);
1275       end
1276     else
1277       % -sp(1) - beta2*Delta_T*[eta_c*psc(1,j,0) + (1/eta_d)*psd(1,j,0)] <= -beta1*(rho*InitialUB+(1-rho)*Initial)/baseMVA
1278       for j = 1:mdi.idx.nj(1)
1279         A = [ -diagBeta2EtaIn1 -diagBeta2overEtaOut1 -Ins ];
1280         u = -beta1(:,1) .* ( rho(:,t).*mdi.Storage.InitialStorageUpperBound + ...
1281                             (1-rho(:,t)).*mdi.Storage.InitialStorage ) / baseMVA;
1282         vs = struct('name', {'Psc', 'Psd', 'Sp'}, 'idx', {{1,j,1}, {1,j,1}, {1}});
1283         om.add_lin_constraint('Sp', {1,j}, A, [], u, vs);
1284       end
1285     end
1286     % Then the rest of the periods
1287     % -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
1288     % where s_I(t,j) = L_I(t,j) * s0 + (Mg(t,j)+Mh(t,j)) * x
1289     for t = 2:nt
1290       for j = 1:mdi.idx.nj(t)
1291         Mj = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :) + ...
1292              mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1293         Lij = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1294         diag1minusRhoBeta1 = spdiags((1-rho(:,t)) .* beta1(:,t), 0, ns, ns);
1295         A = sparse([1:ns,1:ns,1:ns,1:ns]', ...
1296                    [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)]', ...
1297                    [-beta2EtaIn(:,t); -beta2overEtaOut(:,t); beta1(:,t).*rho(:,t); -ones(ns,1)], ...
1298                    ns, nvars) ...
1299                 + diag1minusRhoBeta1 * Mj;
1300         if mdi.Storage.ForceCyclicStorage
1301           As0 = sparse(ns, nvars);
1302           As0(:, vv.i1.S0:vv.iN.S0) = diag1minusRhoBeta1 * Lij;
1303           A = A + As0;
1304           u = zeros(ns, 1);
1305         else
1306           u = full(-diag1minusRhoBeta1 * Lij * mdi.Storage.InitialStorage/baseMVA);
1307         end
1308         om.add_lin_constraint('Sp', {t,j}, A, [], u);
1309       end
1310     end
1311     % Now go on and limit the amount of energy that can be used if a
1312     % contingency does happen. Bound sm first. First examine time period 1 wrt to initial
1313     % stored energy, then t=2 and on.
1314     om.init_indexed_name('lin', 'contSm', {nt, nj_max, nc_max+1});
1315     for j = 1:mdi.idx.nj(1)
1316       for k = 2:mdi.idx.nc(1,j)+1  %% NOTE NO k=1!!!
1317         if mdi.Storage.ForceCyclicStorage
1318           % 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)
1319           A = [ diagBeta4EtaIn1 diagBeta4overEtaOut1 diagBeta3EtaIn1 diagBeta3overEtaOut1 -spdiags(beta5(:,1), 0, ns, ns) ];
1320           u = -MinStorageLevel(:,1)/baseMVA;
1321           vs = struct('name', {'Psc', 'Psd', 'Psc', 'Psd', 'S0'}, ...
1322                       'idx', {{1,j,1}, {1,j,1}, {1,j,k}, {1,j,k}, {}});
1323         else
1324           % 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)
1325           A = [ diagBeta4EtaIn1 diagBeta4overEtaOut1 diagBeta3EtaIn1 diagBeta3overEtaOut1 ];
1326           u = (beta5(:,1).*mdi.Storage.InitialStorageLowerBound - MinStorageLevel(:,1))/baseMVA;
1327           vs = struct('name', {'Psc', 'Psd', 'Psc', 'Psd'}, ...
1328                       'idx', {{1,j,1}, {1,j,1}, {1,j,k}, {1,j,k}});
1329         end
1330         om.add_lin_constraint('contSm', {1,j,k}, A, [], u, vs);
1331       end
1332     end
1333     % then the rest of the periods
1334     % -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)
1335     % where s_I(t,j) = L_I(t,j) * s0 + (Mg(t,j)+Mh(t,j)) * x
1336     for t = 2:nt
1337       for j = 1:mdi.idx.nj(t)
1338         for k = 2:mdi.idx.nc(t,j)+1
1339           Mj = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :) + ...
1340                mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1341           Lij = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1342           diag1minusRhoBeta5 = spdiags((1-rho(:,t)) .* beta5(:,t), 0, ns, ns);
1343           A = sparse([1:ns,1:ns,1:ns,1:ns,1:ns]', ...
1344                      [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)]', ...
1345                      [beta4EtaIn(:,t); beta4overEtaOut(:,t); beta3EtaIn(:,t); beta3overEtaOut(:,t); -beta5(:,t).*rho(:,t)], ...
1346                      ns, nvars) ...
1347                   - diag1minusRhoBeta5 * Mj;
1348           u = -MinStorageLevel(:,t)/baseMVA;
1349           if mdi.Storage.ForceCyclicStorage
1350             As0 = sparse(ns, nvars);
1351             As0(:, vv.i1.S0:vv.iN.S0) = -diag1minusRhoBeta5 * Lij;
1352             A = A + As0;
1353           else
1354             u = u + diag1minusRhoBeta5 * Lij * mdi.Storage.InitialStorageLowerBound/baseMVA;
1355           end
1356           om.add_lin_constraint('contSm', {t,j,k}, A, [], u);
1357         end
1358       end
1359     end
1360     % Bound sp first. First examine time period 1 wrt to initial
1361     % stored energy, then t=2 and on.
1362     om.init_indexed_name('lin', 'contSp', {nt, nj_max, nc_max+1});
1363     for j = 1:mdi.idx.nj(1)
1364       for k = 2:mdi.idx.nc(1,j)+1
1365         if mdi.Storage.ForceCyclicStorage
1366           % -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)
1367           A = [ -diagBeta4EtaIn1 -diagBeta4overEtaOut1 -diagBeta3EtaIn1 -diagBeta3overEtaOut1 spdiags(beta5(:,1), 0, ns, ns)];
1368           u = MaxStorageLevel(:,1)/baseMVA;
1369           vs = struct('name', {'Psc', 'Psd', 'Psc', 'Psd', 'S0'}, ...
1370                       'idx', {{1,j,1}, {1,j,1}, {1,j,k}, {1,j,k}, {}});
1371         else
1372           % -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)
1373           A = [ -diagBeta4EtaIn1 -diagBeta4overEtaOut1 -diagBeta3EtaIn1 -diagBeta3overEtaOut1 ];
1374           u = (MaxStorageLevel(:,1) - beta5(:,1).*mdi.Storage.InitialStorageUpperBound)/baseMVA;
1375           vs = struct('name', {'Psc', 'Psd', 'Psc', 'Psd'}, ...
1376                       'idx', {{1,j,1}, {1,j,1}, {1,j,k}, {1,j,k}});
1377         end
1378         om.add_lin_constraint('contSp', {1,j,k}, A, [], u, vs);
1379       end
1380     end
1381     % then the rest of the periods
1382     % 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)
1383     % where s_I(t,j) = L_I(t,j) * s0 + (Mg(t,j)+Mh(t,j)) * x
1384     for t = 2:nt
1385       for j = 1:mdi.idx.nj(t)
1386         for k = 2:mdi.idx.nc(t,j)+1
1387           Mj = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :) + ...
1388                mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1389           Lij = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1390           diag1minusRhoBeta5 = spdiags((1-rho(:,t)) .* beta5(:,t), 0, ns, ns);
1391           A = sparse([1:ns,1:ns,1:ns,1:ns,1:ns]', ...
1392                      [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)]', ...
1393                      [-beta4EtaIn(:,t); -beta4overEtaOut(:,t); -beta3EtaIn(:,t); -beta3overEtaOut(:,t); beta5(:,t).*rho(:,t)], ...
1394                      ns, nvars) ...
1395                   + diag1minusRhoBeta5 * Mj;
1396           u = MaxStorageLevel(:,t)/baseMVA;
1397           if mdi.Storage.ForceCyclicStorage
1398             As0 = sparse(ns, nvars);
1399             As0(:, vv.i1.S0:vv.iN.S0) = diag1minusRhoBeta5 * Lij;
1400             A = A + As0;
1401           else
1402             u = u - diag1minusRhoBeta5 * Lij * mdi.Storage.InitialStorageUpperBound/baseMVA;
1403           end
1404           om.add_lin_constraint('contSp', {t,j,k}, A, [], u);
1405         end
1406       end
1407     end
1408   end
1409   
1410   % Now, if required, constrain the expected terminal storage quantity; two
1411   % different ways:
1412   if mdi.Storage.ForceExpectedTerminalStorage && mdi.Storage.ForceCyclicStorage
1413     error('most: ForceExpectedTerminalStorage and ForceCyclicStorage cannot be simultaneously true.');
1414   end
1415   if ns
1416     % The following code assumes that no more variables will be added
1417     if mdi.Storage.ForceExpectedTerminalStorage
1418       % 1) Constrain the expected terminal storage to be some target value
1419       A = sparse(ns, nvars);
1420       b = zeros(ns, 1);
1421       for j = 1:mdi.idx.nj(nt)
1422         Ngj = mdi.tstep(nt).Ng( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1423         Nhj = mdi.tstep(nt).Nh( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1424         Lfj = mdi.tstep(nt).Lf( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1425         A = A + mdi.CostWeights(1,j,nt) * (Ngj + Nhj);
1426         b = b + mdi.CostWeights(1,j,nt) * (Lfj * mdi.Storage.InitialStorage) / baseMVA;
1427       end
1428       endprob = sum(mdi.CostWeights(1,1:mdi.idx.nj(nt),nt)');
1429       A = (1/endprob) * A;
1430       b = (1/endprob) * b;
1431       l = mdi.Storage.ExpectedTerminalStorageMin / baseMVA - b;
1432       u = mdi.Storage.ExpectedTerminalStorageMax / baseMVA - b;
1433       om.add_lin_constraint('ESnt', A, l, u);
1434     elseif mdi.Storage.ForceCyclicStorage
1435       % 2) Constrain the initial storage (a variable) to be the same as the final expected storage
1436       A = sparse(ns, nvars);
1437       for j = 1:mdi.idx.nj(nt)
1438         Ngj = mdi.tstep(nt).Ng( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1439         Nhj = mdi.tstep(nt).Nh( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1440         A = A + mdi.CostWeights(1,j,nt) * (Ngj + Nhj);
1441       end
1442       endprob = sum(mdi.CostWeights(1,1:mdi.idx.nj(nt),nt)');
1443       A = (1/endprob) * A;
1444       for j = 1:mdi.idx.nj(nt)
1445         Lfj = mdi.tstep(nt).Lf(j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1446         A(:, vv.i1.S0:vv.iN.S0) = A(:, vv.i1.S0:vv.iN.S0) ...
1447                   + mdi.CostWeights(1,j,nt) * Lfj;
1448       end
1449       A(:, vv.i1.S0:vv.iN.S0) = (1/endprob) * A(:, vv.i1.S0:vv.iN.S0) - speye(ns);
1450       b = zeros(ns, 1);
1451       om.add_lin_constraint('ESnt', A, b, b);
1452     end
1453   end
1454 
1455   % Dynamical system contraints
1456   if nzds || nyds
1457     if verbose
1458       fprintf('  - Building dynamical system constraints.\n');
1459     end
1460     % Compute matrices that give the expected dispatch in time period t
1461     % given that we make it to that period, for all generators at once,
1462     % when multiplied by x, i.e. E[p(t)] = E(t) * x
1463     for t = 1:nt
1464       E = sparse(ng, nvars);
1465       for j = 1:mdi.idx.nj(t)
1466         for k = 1:mdi.idx.nc(t,j)+1;
1467           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);
1468         end
1469       end
1470       mdi.tstep(t).E = E;
1471     end
1472   end
1473 
1474   % Form the dynamical system state equations and bound constraints on the
1475   % state vector
1476   if nzds
1477     om.init_indexed_name('lin', 'DSz', {ntds-1});
1478     b = zeros(nzds, 1);
1479     for t = 1:ntds-1
1480       if t <= nt  % We have p(t) available to drive the dynamical system up to t=nt
1481         % Form the constraint matrix, so B*E*x + A*z(t) - I*z(t+1) = 0
1482         A = mdi.dstep(t).B * mdi.tstep(t).E;
1483       else
1484         % The dynamical system horizon is longer than the injection planning
1485         % horizon and we don't know what p(t) is, but continue to drive the
1486         % dynamical system as if p(t) = 0 and perhaps take that into account
1487         % when setting Ymax, Ymin in this time window.  That is, A*z(t) - I*z(t+1) = 0
1488         A = sparse(nzds, nvars);
1489       end
1490       A(:, vv.i1.Z(t):vv.iN.Z(t)) = mdi.dstep(t).A;
1491       A(:, vv.i1.Z(t+1):vv.iN.Z(t+1)) = -speye(nzds);
1492       om.add_lin_constraint('DSz', {t}, A, b, b);
1493     end
1494   end
1495   
1496   % Form the output equations and their restrictions
1497   if nyds
1498     om.init_indexed_name('lin', 'DSy', {ntds});
1499     for t = 1:ntds
1500       if t <= nt
1501         A = mdi.dstep(t).D * mdi.tstep(t).E;
1502       else
1503         A = sparse(nyds, nvars);
1504       end
1505       if nzds
1506         A(:, vv.i1.Z(t):vv.iN.Z(t)) = mdi.dstep(t).C;
1507       end
1508       l = mdi.dstep(t).ymin;
1509       u = mdi.dstep(t).ymax;
1510       om.add_lin_constraint('DSy', {t}, A, l, u);
1511     end
1512   end
1513   
1514   % UNIT COMMITMENT
1515   if UC
1516     if verbose
1517       fprintf('  - Building unit commitment constraints.\n');
1518     end
1519     % u(t,i) - u(t-1,i) - v(t,i) + w(t,i) = 0
1520     om.init_indexed_name('lin', 'uvw', {nt});
1521     for t = 1:nt
1522       if t == 1
1523         % First for t=1 when u(t-1,i) is really u(0,i) or u(nt,i)
1524         if mdi.UC.CyclicCommitment
1525           vs = struct('name', {'u', 'u', 'v', 'w'}, 'idx', {{1}, {nt}, {1}, {1}});
1526           A = [Ing -Ing -Ing Ing];
1527           b = zeros(ng, 1);
1528         else
1529           vs = struct('name', {'u', 'v', 'w'}, 'idx', {{1}, {1}, {1}});
1530           A = [Ing -Ing Ing];
1531           b = (mdi.UC.InitialState > 0);
1532         end
1533       else
1534         % Then for rest of periods
1535         vs = struct('name', {'u', 'u', 'v', 'w'}, 'idx', {{t}, {t-1}, {t}, {t}});
1536         A = [Ing -Ing -Ing Ing];
1537         b = zeros(ng, 1);
1538       end
1539       om.add_lin_constraint('uvw', {t}, A, b, b, vs);
1540     end
1541     % Then continue with minimimum up time constraints. Again, two
1542     % different forms depending on whether the horizon is cyclical or not
1543     om.init_indexed_name('lin', 'minup', {nt, ng});
1544     for t = 1:nt
1545       for i = 1:ng
1546         ti = t-mdi.UC.MinUp(i)+1:t;
1547         if mdi.UC.CyclicCommitment     % window is circular
1548           for tt = 1:length(ti)
1549             if ti(tt) < 1
1550               ti(tt) = nt + ti(tt);
1551             end
1552           end
1553         end
1554         % limit to positive time
1555         % even with CyclicCommitment, in case MinUp is longer than horizon
1556         % (which implies always ON or always OFF)
1557         ti = ti(ti>0);
1558         vs = struct('name', {'u'}, 'idx', {{t}});
1559         A = sparse(1, i, -1, 1, ng);
1560         for tt = 1:length(ti)
1561             vs(end+1).name = 'v';
1562             vs(end).idx  = {ti(tt)};
1563             A = [A sparse(1, i, 1, 1, ng)];
1564         end
1565         om.add_lin_constraint('minup', {t, i}, A, [], 0, vs);
1566       end
1567     end
1568     % Continue with minimimum downtime constraints. Two
1569     % different forms depending on whether the horizon is cyclical or not
1570     om.init_indexed_name('lin', 'mindown', {nt, ng});
1571     for t = 1:nt
1572       for i = 1:ng
1573         ti = t-mdi.UC.MinDown(i)+1:t;
1574         if mdi.UC.CyclicCommitment     % window is circular
1575           for tt = 1:length(ti)
1576             if ti(tt) < 1
1577               ti(tt) = nt + ti(tt);
1578             end
1579           end
1580         end
1581         % limit to positive time
1582         % even with CyclicCommitment, in case MinDown is longer than horizon
1583         % (which implies always ON or always OFF)
1584         ti = ti(ti>0);
1585         vs = struct('name', {'u'}, 'idx', {{t}});
1586         A = sparse(1, i, 1, 1, ng);
1587         for tt = 1:length(ti)
1588             vs(end+1).name = 'w';
1589             vs(end).idx  = {ti(tt)};
1590             A = [A sparse(1, i, 1, 1, ng)];
1591         end
1592         om.add_lin_constraint('mindown', {t, i}, A, [], 1, vs);
1593       end
1594     end
1595     % Limit generation ranges based on commitment status; first Pmax;
1596     % p - u*Pmax <= 0
1597     % For contingent flows, however, if a generator is ousted as a result
1598     % of the contingency, then this constraint should not be enforced.
1599     om.init_indexed_name('lin', 'uPmax', {nt, nj_max, nc_max+1});
1600     for t = 1:nt
1601       for j = 1:mdi.idx.nj(t)
1602         for k = 1:mdi.idx.nc(t,j)+1
1603           mpc = mdi.flow(t,j,k).mpc;
1604           ii = find(mpc.gen(:, GEN_STATUS));
1605           nii = length(ii);
1606           vs = struct('name', {'Pg', 'u'}, 'idx', {{t,j,k}, {t}});
1607           A = [ sparse(1:nii, ii, 1, nii, ng) ...
1608                 sparse(1:nii, ii, -mpc.gen(ii, PMAX)/baseMVA, nii, ng) ];
1609           u = zeros(nii, 1);
1610           om.add_lin_constraint('uPmax', {t,j,k}, A, [], u, vs);
1611         end
1612       end
1613     end
1614     % Then Pmin,  -p + u*Pmin <= 0
1615     om.init_indexed_name('lin', 'uPmin', {nt, nj_max, nc_max+1});
1616     for t = 1:nt
1617       for j = 1:mdi.idx.nj(t)
1618         for k = 1:mdi.idx.nc(t,j)+1
1619           mpc = mdi.flow(t,j,k).mpc;
1620           ii = find(mpc.gen(:, GEN_STATUS));
1621           nii = length(ii);
1622           vs = struct('name', {'Pg', 'u'}, 'idx', {{t,j,k}, {t}});
1623           A = [ sparse(1:nii, ii, -1, nii, ng) ...
1624                 sparse(1:nii, ii, mpc.gen(ii, PMIN)/baseMVA, nii, ng) ];
1625           u = zeros(nii, 1);
1626           om.add_lin_constraint('uPmin', {t,j,k}, A, [], u, vs);
1627         end
1628       end
1629     end
1630     % Then, if there is Qg coordination, do the same for Qg
1631     % q - u*Qmax <= 0
1632     % For contingent flows, however, if a generator is ousted as a result
1633     % of the contingency, then this constraint should not be enforced.
1634     if mdi.QCoordination
1635       om.init_indexed_name('lin', 'uQmax', {nt, nj_max, nc_max+1});
1636       for t = 1:nt
1637         for j = 1:mdi.idx.nj(t)
1638           for k = 1:mdi.idx.nc(t,j)+1
1639             mpc = mdi.flow(t,j,k).mpc;
1640             ii = find(mpc.gen(:, GEN_STATUS));
1641             nii = length(ii);
1642             vs = struct('name', {'Qg', 'u'}, 'idx', {{t,j,k}, {t}});
1643             A = [ sparse(1:nii, ii, 1, nii, ng) ...
1644                   sparse(1:nii, ii, -mpc.gen(ii, QMAX)/baseMVA, nii, ng) ];
1645             u = zeros(nii, 1);
1646             om.add_lin_constraint('uQmax', {t,j,k}, A, [], u, vs);
1647           end
1648         end
1649       end
1650       % Then Qmin,  -q + u*Qmin <= 0
1651       om.init_indexed_name('lin', 'uQmin', {nt, nj_max, nc_max+1});
1652       for t = 1:nt
1653         for j = 1:mdi.idx.nj(t)
1654           for k = 1:mdi.idx.nc(t,j)+1
1655             mpc = mdi.flow(t,j,k).mpc;
1656             ii = find(mpc.gen(:, GEN_STATUS));
1657             nii = length(ii);
1658             vs = struct('name', {'Qg', 'u'}, 'idx', {{t,j,k}, {t}});
1659             A = [ sparse(1:nii, ii, -1, nii, ng) ...
1660                   sparse(1:nii, ii, mpc.gen(ii, QMIN)/baseMVA, nii, ng) ];
1661             u = zeros(nii, 1);
1662             om.add_lin_constraint('uQmin', {t,j,k}, A, [], u, vs);
1663           end
1664         end
1665       end
1666     end
1667   end
1668 
1669   if verbose
1670     fprintf('- Building cost structures.\n');
1671   end
1672   % Start building the cost.  Two main components, the input data cost and
1673   % the coordination cost are specified.  The coordination cost is assumed to
1674   % have been buit with knowledge of the variable structure, and is simply
1675   % passed on.  The input data cost is assembled into the appropriate
1676   % spots.
1677   %
1678   % f = 0.5 * x' * (H1 + Hcoord) * x + (C1' + Ccoord) * x + c1 + ccoord
1679 
1680   % First assign the ramping costs; H1 has few coefficients initially and
1681   % this should make the shuffling and reordering of coefficients more
1682   % efficient.  All other accesses to H1 will be diagonal insertions, which
1683   % take less time than anti-diagonal insertions.
1684   % First do first period wrt to InitialPg.
1685   if mdi.OpenEnded
1686     om.init_indexed_name('qdc', 'RampWear', {nt, nj_max, nj_max});
1687   else
1688     om.init_indexed_name('qdc', 'RampWear', {nt+1, nj_max, nj_max});
1689   end
1690   for j = 1:mdi.idx.nj(1)
1691     w = mdi.tstep(1).TransMat(j,1);  % the probability of going from initial state to jth
1692     Q = spdiags(w * baseMVA^2 * mdi.RampWearCostCoeff(:,1), 0, ng, ng);
1693     c = -w * baseMVA * mdi.RampWearCostCoeff(:,1) .* mdi.InitialPg;
1694     vs = struct('name', {'Pg'}, 'idx', {{1,j,1}});
1695     k0 = w * 0.5 * mdi.RampWearCostCoeff(:,1)' * mdi.InitialPg.^2;
1696     om.add_quad_cost('RampWear', {1,j,1}, Q, c, k0, vs);
1697   end
1698   % Then the remaining periods
1699   for t = 2:nt
1700     for j2 = 1:mdi.idx.nj(t)
1701       for j1 = 1:mdi.idx.nj(t-1)
1702         w = mdi.tstep(t).TransMat(j2,j1) * mdi.CostWeights(1, j1, t-1);
1703         h = w * baseMVA^2 * mdi.RampWearCostCoeff(:,t);
1704         i = (1:ng)';
1705         j = ng+(1:ng)';
1706         Q = sparse([i;j;i;j], [i;i;j;j], [h;-h;-h;h], 2*ng, 2*ng);
1707         vs = struct('name', {'Pg', 'Pg'}, 'idx', {{t-1,j1,1}, {t,j2,1}});
1708         om.add_quad_cost('RampWear', {t,j1,j2}, Q, zeros(2*ng,1), 0, vs);
1709       end
1710     end
1711   end
1712   % Finally, if there is a terminal state problem, apply cost to
1713   % the transition starting from t=nt.  Note that in this case
1714   % mdi.tstep(nt+1).TransMat must be defined! it is the only piece of data
1715   % that makes sense for nt+1; all other fields in mdi.tstep(nt+1) can be empty.
1716   if ~mdi.OpenEnded
1717     for j = 1:mdi.idx.nj(nt)
1718       w = mdi.tstep(nt+1).TransMat(1, j) * mdi.CostWeights(1, j, nt);
1719       Q = spdiags(w * baseMVA^2 * mdi.RampWearCostCoeff(:,nt+1), 0, ng, ng);
1720       c = -w * baseMVA * mdi.RampWearCostCoeff(:,nt+1) .* mdi.TerminalPg;
1721       vs = struct('name', {'Pg'}, 'idx', {{nt,j,1}});
1722       k0 = w * 0.5 * mdi.RampWearCostCoeff(:,nt+1)' * mdi.TerminalPg.^2;
1723       om.add_quad_cost('RampWear', {nt+1,j,1}, Q, c, k0, vs);
1724     end
1725   end
1726 
1727   % Now go on and assign energy, inc/dec and contingency reserves
1728   % costs for all committed units.
1729   om.init_indexed_name('qdc', 'Cp', {nt, nj_max, nc_max+1});
1730   om.init_indexed_name('qdc', 'Cy', {nt, nj_max, nc_max+1});
1731   om.init_indexed_name('qdc', 'Cpp', {nt, nj_max, nc_max+1});
1732   om.init_indexed_name('qdc', 'Cpm', {nt, nj_max, nc_max+1});
1733   if mdi.IncludeFixedReserves
1734     om.init_indexed_name('qdc', 'Rcost', {nt, nj_max, nc_max+1});
1735   end
1736   om.init_indexed_name('qdc', 'Crpp', {nt});
1737   om.init_indexed_name('qdc', 'Crpm', {nt});
1738   for t = 1:nt
1739     for j = 1:mdi.idx.nj(t)
1740       for k = 1:mdi.idx.nc(t,j)+1
1741         w = mdi.CostWeightsAdj(k,j,t);     %% NOTE (k,j,t) order !!!
1742 
1743         % weighted polynomial energy costs for committed units
1744         gc = mdi.flow(t,j,k).mpc.gencost;
1745         ipol = find(gc(:, MODEL) == POLYNOMIAL);
1746         if ~isempty(ipol)
1747           ncost = gc(ipol(1), NCOST);
1748           if all(gc(ipol, NCOST) == ncost)    %% uniform order of polynomials
1749             %% use vectorized code
1750             if ncost > 3
1751               error('most: polynomial generator costs of order higher than quadratic not supported');
1752             elseif ncost == 3
1753               Q = sparse(ipol, ipol, 2 * w * baseMVA^2*gc(ipol, COST), ng, ng);
1754             else
1755               Q = sparse(ng,ng);
1756             end
1757             c = zeros(ng, 1);
1758             if ncost >= 2
1759               c(ipol) = w * baseMVA*gc(ipol, COST+ncost-2);
1760             end
1761             k0 = w * sum(gc(ipol, COST+ncost-1));
1762           else                                %% non-uniform order of polynomials
1763             %% use a loop
1764             Q = sparse(ng,ng);
1765             c = zeros(ng, 1);
1766             for i = ipol'
1767               ncost = gc(i, NCOST);
1768               if ncost > 3
1769                 error('most: polynomial generator costs of order higher than quadratic not supported');
1770               elseif ncost == 3
1771                 Q(i,i) = 2 * w * baseMVA^2*gc(i, COST);
1772               end
1773               if ncost >= 2
1774                 c(i) = w * baseMVA*gc(i, COST+ncost-2);
1775               end
1776               k0 = w * gc(i, COST+ncost-1);
1777             end
1778           end
1779           vs = struct('name', {'Pg'}, 'idx', {{t,j,k}});
1780           om.add_quad_cost('Cp', {t,j,k}, Q, c, k0, vs);
1781         end
1782 
1783         % weighted y-variables for piecewise linear energy costs for committed units
1784         % ipwl = find( (mdi.flow(t,j,k).mpc.gen(:,GEN_STATUS) > 0) & (gc(:,MODEL) == PW_LINEAR));
1785         if mdi.idx.ny(t,j,k)
1786           c = w * ones(mdi.idx.ny(t,j,k),1);
1787           vs = struct('name', {'y'}, 'idx', {{t,j,k}});
1788           om.add_quad_cost('Cy', {t,j,k}, [], c, 0, vs);
1789         end
1790 
1791         % inc and dec offers for each flow
1792         c = w * baseMVA * mdi.offer(t).PositiveActiveDeltaPrice(:);
1793         vs = struct('name', {'dPp'}, 'idx', {{t,j,k}});
1794         om.add_quad_cost('Cpp', {t,j,k}, [], c, 0, vs);
1795         c = w * baseMVA * mdi.offer(t).NegativeActiveDeltaPrice(:);
1796         vs = struct('name', {'dPm'}, 'idx', {{t,j,k}});
1797         om.add_quad_cost('Cpm', {t,j,k}, [], c, 0, vs);
1798 
1799         % weighted fixed reserves cost
1800         if mdi.IncludeFixedReserves
1801           c = w * mdi.FixedReserves(t,j,k).cost(r.igr) * baseMVA;
1802           vs = struct('name', {'R'}, 'idx', {{t,j,k}});
1803           om.add_quad_cost('Rcost', {t,j,k}, [], c, 0, vs);
1804         end
1805       end
1806     end
1807     
1808     % contingency reserve costs
1809     c = baseMVA * mdi.StepProb(t) * mdi.offer(t).PositiveActiveReservePrice(:);
1810     vs = struct('name', {'Rpp'}, 'idx', {{t}});
1811     om.add_quad_cost('Crpp', {t}, [], c, 0, vs);
1812     c = baseMVA * mdi.StepProb(t) * mdi.offer(t).NegativeActiveReservePrice(:);
1813     vs = struct('name', {'Rpm'}, 'idx', {{t}});
1814     om.add_quad_cost('Crpm', {t}, [], c, 0, vs);
1815   end
1816   % Assign load following ramp reserve costs.  Do first nt-1 periods first
1817   om.init_indexed_name('qdc', 'Crrp', {mdi.idx.ntramp});
1818   om.init_indexed_name('qdc', 'Crrm', {mdi.idx.ntramp});
1819   for t = 1:nt-1,
1820     c = baseMVA * mdi.StepProb(t+1) * mdi.offer(t).PositiveLoadFollowReservePrice(:);
1821     vs = struct('name', {'Rrp'}, 'idx', {{t}});
1822     om.add_quad_cost('Crrp', {t}, [], c, 0, vs);
1823     c = baseMVA * mdi.StepProb(t+1) * mdi.offer(t).NegativeLoadFollowReservePrice(:);
1824     vs = struct('name', {'Rrm'}, 'idx', {{t}});
1825     om.add_quad_cost('Crrm', {t}, [], c, 0, vs);
1826   end
1827   % Then do last period if needed Terminal state case
1828   if ~mdi.OpenEnded
1829     %% are these costs missing a mdi.StepProb(t)?  -- rdz
1830     c = baseMVA * mdi.offer(nt).PositiveLoadFollowReservePrice(:);
1831     vs = struct('name', {'Rrp'}, 'idx', {{nt}});
1832     om.add_quad_cost('Crrp', {nt}, [], c, 0, vs);
1833     c = baseMVA * mdi.offer(nt).NegativeLoadFollowReservePrice(:);
1834     vs = struct('name', {'Rrm'}, 'idx', {{nt}});
1835     om.add_quad_cost('Crrm', {nt}, [], c, 0, vs);
1836   end
1837   % Assign startup/shutdown costs, if any, and fixed operating costs
1838   if UC
1839     om.init_indexed_name('qdc', 'c00', {nt});
1840     om.init_indexed_name('qdc', 'startup', {nt});
1841     om.init_indexed_name('qdc', 'shutdown', {nt});
1842     for t = 1:nt
1843       ww = zeros(ng, 1);
1844       for j = 1:mdi.idx.nj(t)
1845         for k = 1:mdi.idx.nc(t)+1
1846           ww = ww + mdi.CostWeightsAdj(k,j,t) * mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS);
1847         end
1848       end
1849       c = ww.*mdi.UC.c00(:,t);
1850       vs = struct('name', {'u'}, 'idx', {{t}});
1851       om.add_quad_cost('c00', {t}, [], c, 0, vs);
1852       c = mdi.StepProb(t)*mdi.flow(t,1,1).mpc.gencost(:, STARTUP);
1853       vs = struct('name', {'v'}, 'idx', {{t}});
1854       om.add_quad_cost('startup', {t}, [], c, 0, vs);
1855       c = mdi.StepProb(t)*mdi.flow(t,1,1).mpc.gencost(:, SHUTDOWN);
1856       vs = struct('name', {'w'}, 'idx', {{t}});
1857       om.add_quad_cost('shutdown', {t}, [], c, 0, vs);
1858     end
1859   end
1860   % Finally, assign any value to leftover stored energy
1861   if ns
1862     A1 = sparse(ns, ns);
1863     A2 = sparse(ns, nvars);
1864     A3 = sparse(ns, nvars);
1865     A4 = sparse(ns, nvars);
1866     A5 = sparse(ns, nvars);
1867     A6 = sparse(ns, nvars);
1868     A7 = sparse(ns, nvars);
1869     % The following code assumes that no more variables will be added
1870     vv = om.get_idx();
1871     for t = 1:nt
1872       % Compute cost coefficients for value of expected leftover storage
1873       % after a contingency
1874       for j = 1:mdi.idx.nj(t)
1875         % pick rows for jth base injections
1876         Gtj0  = mdi.tstep(t).G(  j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1877         Htj0  = mdi.tstep(t).H(  j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1878         Litj0 = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1879         Mgtj0 = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1880         Mhtj0 = mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1881         sum_psi_tjk = sum(mdi.CostWeights(2:mdi.idx.nc(t,j)+1,j,t));
1882         if t == nt
1883           A1 = A1 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta1(:,t), 0, ns, ns) * Litj0;
1884           A2 = A2 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta1(:,t), 0, ns, ns) * Mgtj0;
1885           A3 = A3 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta1(:,t), 0, ns, ns) * Mhtj0;
1886           A4 = A4 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta2(:,t), 0, ns, ns) * Gtj0;
1887           A5 = A5 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta2(:,t), 0, ns, ns) * Htj0;
1888         end
1889         A1 = A1 + sum_psi_tjk * spdiags(OutEff(:,t) .* beta5(:,t), 0, ns, ns) * Litj0;
1890         A2 = A2 + sum_psi_tjk * (spdiags(OutEff(:,t) .* beta5(:,t), 0, ns, ns) * Mgtj0 + spdiags(OutEff(:,t) .* beta4(:,t), 0, ns, ns) * Gtj0);
1891         A3 = A3 + sum_psi_tjk * (spdiags(OutEff(:,t) .* beta5(:,t), 0, ns, ns) * Mhtj0 + spdiags(OutEff(:,t) .* beta4(:,t), 0, ns, ns) * Htj0);
1892         for k = 2:mdi.idx.nc(t,j)+1
1893           ii  = (1:ns)';
1894           jj1 = (vv.i1.Psc(t,j,k):vv.iN.Psc(t,j,k))';
1895           jj2 = (vv.i1.Psd(t,j,k):vv.iN.Psd(t,j,k))';
1896           Gtjk = sparse(ii, jj1, -mdi.Delta_T  *  InEff(:,t), ns, nvars);
1897           Htjk = sparse(ii, jj2, -mdi.Delta_T ./ OutEff(:,t), ns, nvars);
1898           A6 = A6 + mdi.CostWeights(k,j,t) * spdiags(OutEff(:,t) .* beta3(:,t), 0, ns, ns) * Gtjk;
1899           A7 = A7 + mdi.CostWeights(k,j,t) * spdiags(OutEff(:,t) .* beta3(:,t), 0, ns, ns) * Htjk;
1900         end
1901       end
1902     end
1903     Cfstor = -baseMVA * ...
1904        (mdi.Storage.TerminalStoragePrice'      * (A2 + A3) + ...
1905         mdi.Storage.TerminalChargingPrice0'    * A4 + ...
1906         mdi.Storage.TerminalDischargingPrice0' * A5 + ...
1907         mdi.Storage.TerminalChargingPriceK'    * A6 + ...
1908         mdi.Storage.TerminalDischargingPriceK' * A7);
1909     if mdi.Storage.ForceCyclicStorage
1910       % If the horizon model for the storage is cyclic and therefore s0 is a
1911       % variable, then that initial storage must come at a cost,
1912       % (InitialStorageCost) otherwise the optimizer will force the gratis
1913       % s0 up just to have (possibly) more storage left at the end.
1914       Cfstor(vv.i1.S0:vv.iN.S0) = ...
1915         Cfstor(vv.i1.S0:vv.iN.S0) + ...
1916             baseMVA * mdi.Storage.InitialStorageCost';
1917       % and the term in the final expected storage related to s0 is also
1918       % not constant, so must be included in the objective function
1919       Cfstor(vv.i1.S0:vv.iN.S0) = ...
1920           Cfstor(vv.i1.S0:vv.iN.S0) - ...
1921           baseMVA * mdi.Storage.TerminalStoragePrice' * A1;
1922     end
1923     om.add_quad_cost('fstor', [], Cfstor', 0);
1924 
1925     % The following is a hack to make the storage state bounds tight;
1926     % assign them a very small cost
1927     om.init_indexed_name('qdc', 'SpSmFudge', {nt});
1928     c = 1e-2 * [-ones(ns,1); ones(ns,1)];
1929     for t = 1:nt
1930       vs = struct('name', {'Sm', 'Sp'}, 'idx', {{t}, {t}});
1931       om.add_quad_cost('SpSmFudge', {t}, [], c, 0, vs);
1932     end
1933   else
1934     Cfstor = sparse(1, nvars);
1935   end
1936 
1937   %% handing of user-defined constraints and costs would go here
1938 
1939   % Asssemble contraints, variable bounds and costs
1940   if verbose
1941     fprintf('- Assembling full set of constraints.\n');
1942   end
1943   [mdi.QP.A, mdi.QP.l, mdi.QP.u] = om.params_lin_constraint();
1944   if verbose
1945     fprintf('- Assembling full set of variable bounds.\n');
1946   end
1947   [mdi.QP.x0, mdi.QP.xmin, mdi.QP.xmax, mdi.QP.vtype] = om.params_var();
1948   if verbose
1949     fprintf('- Assembling full set of costs.\n');
1950   end
1951   [mdi.QP.H1, mdi.QP.C1, mdi.QP.c1] = om.params_quad_cost();
1952 
1953   % Plug into struct
1954   mdi.QP.Cfstor = Cfstor;
1955   mdi.om = om;
1956 else
1957   om = mdi.om;
1958 end     % if mpopt.most.build_model
1959 
1960 % With all pieces of the cost in place, can proceed to build the total
1961 % cost now.
1962 if isfield(mdi, 'CoordCost') && ...
1963         (~isempty(mdi.CoordCost.Cuser) || ~isempty(mdi.CoordCost.Huser))
1964   if verbose
1965     fprintf('- Adding coordination cost to standard cost.\n');
1966   end
1967   nvuser = length(mdi.CoordCost.Cuser);
1968   nvdiff = mdi.idx.nvars - nvuser;
1969   om.add_quad_cost( 'CoordCost', ...
1970                     [ mdi.CoordCost.Huser    sparse(nvuser,nvdiff);
1971                       sparse(nvdiff,nvuser)  sparse(nvdiff,nvdiff) ], ...
1972                     mdi.CoordCost.Cuser(:), mdi.CoordCost.cuser);
1973 end
1974 [mdi.QP.H, mdi.QP.C, mdi.QP.c] = om.params_quad_cost();
1975 
1976 et_setup = toc(t0);
1977 t0 = tic;
1978 
1979 % Call solver!
1980 mdo = mdi;              %% initialize output
1981 mdo.om = om.copy();     %% make copy of opt_model object, so changes to
1982                         %% output obj (mdo) don't modify input obj (mdi)
1983 if mpopt.most.solve_model
1984   %% check consistency of model options (in case mdi was built in previous call)
1985   if mdi.DCMODEL ~= mo.DCMODEL
1986     error('MDI.DCMODEL inconsistent with MPOPT.most.dc_model');
1987   end
1988   if mdi.IncludeFixedReserves ~= mo.IncludeFixedReserves
1989     error('MDI.IncludeFixedReserves inconsistent with MPOPT.most.fixed_res (and possible presence of MDI.FixedReserves(t,j,k))');
1990   end
1991   if mdi.SecurityConstrained ~= mo.SecurityConstrained
1992     error('MDI.SecurityConstrained inconsistent with MPOPT.most.security_constraints (and possible presence of MDI.cont(t,j).contab)');
1993   end
1994   if mdi.QCoordination ~= mo.QCoordination
1995     error('MDI.QCoordination inconsistent with MPOPT.most.q_coordination');
1996   end
1997   if mdi.Storage.ForceCyclicStorage ~= mo.ForceCyclicStorage
1998     error('MDI.Storage.ForceCyclicStorage inconsistent with MPOPT.most.storage.cyclic');
1999   end
2000   if mdi.Storage.ForceExpectedTerminalStorage ~= mo.ForceExpectedTerminalStorage
2001     error('MDI.Storage.ForceExpectedTerminalStorage inconsistent with MPOPT.most.storage.terminal_target (and possible presence of MDI.Storage.ExpectedTerminalStorageAim|Min|Max)');
2002   end
2003   if mdi.UC.run ~= UC
2004     error('MDI.UC.run inconsistent with MPOPT.most.uc.run (and possible presence of MDI.UC.CommitKey)');
2005   end
2006   %% set options
2007   model = om.problem_type();
2008   mdo.QP.opt = mpopt2qpopt(mpopt, model, 'most');
2009   mdo.QP.opt.x0 = [];
2010   if verbose
2011     fprintf('- Calling %s solver.\n\n', model);
2012     fprintf('============================================================================\n\n');
2013   end
2014   [mdo.QP.x, mdo.QP.f, mdo.QP.exitflag, mdo.QP.output, mdo.QP.lambda ] = ...
2015         mdo.om.solve(mdo.QP.opt);
2016   if mdo.QP.exitflag > 0
2017     success = 1;
2018     if verbose
2019       fprintf('\n============================================================================\n');
2020       fprintf('- MOST: %s solved successfully.\n', model);
2021     end
2022   else
2023     success = 0;
2024     if verbose
2025       fprintf('\n============================================================================\n');
2026       fprintf('- MOST: %s solver ''%s'' failed with exit flag = %d\n', model, mdo.QP.opt.alg, mdo.QP.exitflag);
2027     end
2028 %     fprintf('\n============================================================================\n');
2029 %     fprintf('- MOST: %s solver ''%s'' failed with exit flag = %d\n', model, mdo.QP.opt.alg, mdo.QP.exitflag);
2030 %     fprintf('  You can query the workspace to debug.\n')
2031 %     fprintf('  When finished, type the word "dbcont" to continue.\n\n');
2032 %     keyboard;
2033   end
2034   % Unpack results
2035   if verbose
2036     fprintf('- Post-processing results.\n');
2037   end
2038   if success
2039     om = mdo.om;
2040     for t = 1:nt
2041       if UC
2042         mdo.UC.CommitSched(:, t) = om.get_soln('var', 'u', {t});
2043       end
2044       for j = 1:mdi.idx.nj(t)
2045         for k = 1:mdi.idx.nc(t,j)+1
2046           mpc = mdo.flow(t,j,k).mpc;      %% pull mpc from output struct
2047           % Some initialization of data
2048           if mdo.DCMODEL
2049             mpc.bus(:, VM) = 1;
2050           end
2051           % Injections and shadow prices
2052           mpc.gen(:, PG) = baseMVA * om.get_soln('var', 'Pg', {t,j,k});
2053           %% need to update Qg for loads consistent w/constant power factor
2054           Pmin = mpc.gen(:, PMIN);
2055           Qmin = mpc.gen(:, QMIN);
2056           Qmax = mpc.gen(:, QMAX);
2057           ivl = find( isload(mpc.gen) & (Qmin ~= 0 | Qmax ~= 0) );
2058           Qlim = (Qmin(ivl) == 0) .* Qmax(ivl) + (Qmax(ivl) == 0) .* Qmin(ivl);
2059           mpc.gen(ivl, QG) = mpc.gen(ivl, PG) .* Qlim ./ Pmin(ivl);
2060           if mdo.DCMODEL
2061             %% bus angles
2062             Va = om.get_soln('var', 'Va', {t,j,k});
2063             mpc.bus(:, VA) = (180/pi) * Va;
2064           
2065             %% nodal prices
2066             price = (om.get_soln('lin', 'mu_u', 'Pmis', {t,j,k})-om.get_soln('lin', 'mu_l', 'Pmis', {t,j,k})) / baseMVA;
2067             mpc.bus(:, LAM_P) = price;
2068           
2069             %% line flows and line limit shadow prices
2070             mpc.branch(:, PF) = 0;
2071             mpc.branch(:, QF) = 0;
2072             mpc.branch(:, PT) = 0;
2073             mpc.branch(:, QT) = 0;
2074             mpc.branch(:, MU_SF) = 0;
2075             mpc.branch(:, MU_ST) = 0;
2076             ion = find(mpc.branch(:, BR_STATUS));
2077             AA = om.params_lin_constraint('Pf', {t,j,k});
2078             lf = baseMVA * (AA * Va + mdo.flow(t,j,k).PLsh);
2079             mpc.branch(ion, PF) = lf;
2080             mpc.branch(ion, PT) = -lf;
2081             mpc.branch(ion, MU_SF) = om.get_soln('lin', 'mu_u', 'Pf', {t,j,k}) / baseMVA;
2082             mpc.branch(ion, MU_ST) = om.get_soln('lin', 'mu_l', 'Pf', {t,j,k}) / baseMVA;
2083           else
2084             %% system price
2085             price = (om.get_soln('lin', 'mu_l', 'Pmis', {t,j,k})-om.get_soln('lin', 'mu_u', 'Pmis', {t,j,k})) / baseMVA;
2086             mpc.bus(:, LAM_P) = price;
2087           end
2088           if UC
2089             % igenon does not contain gens ousted because of a contingency or
2090             % a forced-off UC.CommitKey
2091             igenon = find(mpc.gen(:, GEN_STATUS));
2092             mpc.gen(igenon, GEN_STATUS) = mdo.UC.CommitSched(igenon, t);
2093             gs = mpc.gen(igenon, GEN_STATUS) > 0; % gen status
2094             mpc.gen(:, MU_PMAX) = 0;
2095             mpc.gen(:, MU_PMIN) = 0;
2096             mpc.gen(igenon, MU_PMAX) = gs .* ...
2097                     om.get_soln('lin', 'mu_u', 'uPmax', {t,j,k}) / baseMVA;
2098             mpc.gen(igenon, MU_PMIN) = gs .* ...
2099                     om.get_soln('lin', 'mu_u', 'uPmin', {t,j,k}) / baseMVA;
2100             if mdo.QCoordination
2101               mpc.gen(:, MU_QMAX) = 0;
2102               mpc.gen(:, MU_QMIN) = 0;
2103               mpc.gen(igenon, MU_QMAX) = gs .* ...
2104                       om.get_soln('lin', 'mu_u', 'uQmax', {t,j,k}) / baseMVA;
2105               mpc.gen(igenon, MU_QMIN) = gs .* ...
2106                       om.get_soln('lin', 'mu_u', 'uQmin', {t,j,k}) / baseMVA;
2107             end
2108           else
2109             gs = mpc.gen(:, GEN_STATUS) > 0;      % gen status
2110             mpc.gen(:, MU_PMAX) = gs .* ...
2111                     om.get_soln('var', 'mu_u', 'Pg', {t,j,k}) / baseMVA;
2112             mpc.gen(:, MU_PMIN) = gs .* ...
2113                     om.get_soln('var', 'mu_l', 'Pg', {t,j,k}) / baseMVA;
2114             if mdo.QCoordination
2115               mpc.gen(:, MU_QMAX) = gs .* ...
2116                       om.get_soln('var', 'mu_u', 'Qg', {t,j,k}) / baseMVA;
2117               mpc.gen(:, MU_QMIN) = gs .* ...
2118                       om.get_soln('var', 'mu_l', 'Qg', {t,j,k}) / baseMVA;
2119             end
2120           end
2121           if mdi.IncludeFixedReserves
2122             z = zeros(ng, 1);
2123             r = mdo.FixedReserves(t,j,k);
2124             r.R   = z;
2125             r.prc = z;
2126             r.mu = struct('l', z, 'u', z, 'Pmax', z);
2127             r.totalcost = sum(om.get_soln('qdc', 'Rcost', {t,j,k}));
2128             r.R(r.igr) = om.get_soln('var', 'R', {t,j,k}) * baseMVA;
2129             R_mu_l = om.get_soln('lin', 'mu_l', 'Rreq', {t,j,k});
2130             for gg = r.igr
2131               iz = find(r.zones(:, gg));
2132               r.prc(gg) = sum(R_mu_l(iz)) / baseMVA;
2133             end
2134             r.mu.l(r.igr) = om.get_soln('var', 'mu_l', 'R', {t,j,k}) / baseMVA;
2135             r.mu.u(r.igr) = om.get_soln('var', 'mu_u', 'R', {t,j,k}) / baseMVA;
2136             r.mu.Pmax(r.igr) = om.get_soln('lin', 'mu_u', 'Pg_plus_R', {t,j,k}) / baseMVA;
2137             mpc.reserves = r;
2138           end
2139           mdo.flow(t,j,k).mpc = mpc;     %% stash modified mpc in output struct
2140         end
2141       end
2142       % Contract, contingency reserves, energy limits
2143       mdo.results.Pc(:,t)  = baseMVA * om.get_soln('var', 'Pc', {t});
2144       mdo.results.Rpp(:,t) = baseMVA * om.get_soln('var', 'Rpp', {t});
2145       mdo.results.Rpm(:,t) = baseMVA * om.get_soln('var', 'Rpm', {t});
2146       if ns
2147         mdo.results.Sm(:,t)  = baseMVA * om.get_soln('var', 'Sm', {t});
2148         mdo.results.Sp(:,t)  = baseMVA * om.get_soln('var', 'Sp', {t});
2149       end
2150     end
2151     % Ramping reserves
2152     for t = 1:mdo.idx.ntramp
2153       mdo.results.Rrp(:,t) = baseMVA * om.get_soln('var', 'Rrp', {t});
2154       mdo.results.Rrm(:,t) = baseMVA * om.get_soln('var', 'Rrm', {t});
2155     end
2156     % Expected energy prices for generators, per generator and per period,
2157     % both absolute and conditional on making it to that period
2158     mdo.results.GenPrices = zeros(ng, nt);
2159     mdo.results.CondGenPrices = zeros(ng, nt);
2160     for t = 1:nt
2161       pp = zeros(ng,1);
2162       for j = 1:mdo.idx.nj(t)
2163         for k = 1:mdo.idx.nc(t,j)+1
2164           pp = pp + mdo.flow(t,j,k).mpc.bus(mdo.flow(t,j,k).mpc.gen(:,GEN_BUS), LAM_P);
2165         end
2166       end
2167       mdo.results.GenPrices(:,t) = pp;
2168       mdo.results.CondGenPrices(:, t) = pp / mdo.StepProb(t);
2169     end
2170     % Obtain contingency reserve prices, per generator and period
2171     mdo.results.RppPrices = zeros(ng, nt);
2172     mdo.results.RpmPrices = zeros(ng, nt);
2173     for t = 1:nt
2174       mdo.results.RppPrices(:, t) = om.get_soln('var', 'mu_l', 'Rpp', {t}) / baseMVA;
2175       mdo.results.RpmPrices(:, t) = om.get_soln('var', 'mu_l', 'Rpm', {t}) / baseMVA;
2176       for j = 1:mdi.idx.nj(t);
2177         for k = 1:mdi.idx.nc(t,j)+1
2178           ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
2179           mdo.results.RppPrices(ii, t) = mdo.results.RppPrices(ii, t) + om.get_soln('lin', 'mu_l', 'dPpRp', {t,j,k}) / baseMVA;
2180           mdo.results.RpmPrices(ii, t) = mdo.results.RpmPrices(ii, t) + om.get_soln('lin', 'mu_l', 'dPmRm', {t,j,k}) / baseMVA;
2181         end
2182       end
2183     end
2184     % Obtain ramping reserve prices, per generator and period
2185     mdo.results.RrpPrices = zeros(ng, mdo.idx.ntramp);
2186     mdo.results.RrmPrices = zeros(ng, mdo.idx.ntramp);
2187     % First, 1:nt-1
2188     for t = 1:nt-1
2189       for j1 = 1:mdo.idx.nj(t)
2190         for j2 = 1:mdo.idx.nj(t+1)
2191           if mdi.tstep(t+1).TransMask(j2,j1)
2192             mdo.results.RrpPrices(:, t) = mdo.results.RrpPrices(:, t) + om.get_soln('lin', 'mu_l', 'Rrp', {t,j1,j2}) / baseMVA;
2193             mdo.results.RrmPrices(:, t) = mdo.results.RrmPrices(:, t) + om.get_soln('lin', 'mu_l', 'Rrm', {t,j1,j2}) / baseMVA;
2194           end
2195         end
2196       end
2197     end
2198     % then last period only if specified for with terminal state
2199     if ~mdo.OpenEnded
2200       for j1 = 1:mdo.idx.nj(nt)
2201         mdo.results.RrpPrices(:, nt) = mdo.results.RrpPrices(:, nt) + om.get_soln('lin', 'mu_l', 'Rrp', {nt,j1,1}) / baseMVA;
2202         mdo.results.RrmPrices(:, nt) = mdo.results.RrmPrices(:, nt) + om.get_soln('lin', 'mu_l', 'Rrm', {nt,j1,1}) / baseMVA;
2203       end
2204     end
2205     % Expected wear and tear costs per gen and period
2206     mdo.results.ExpectedRampCost = zeros(ng, mdo.idx.ntramp+1);
2207     % First do first period wrt to InitialPg.
2208     for j = 1:mdi.idx.nj(1)
2209       w = mdo.tstep(1).TransMat(j,1); % the probability of going from initial state to jth
2210       mdo.results.ExpectedRampCost(:, 1) = mdo.results.ExpectedRampCost(:, 1) ...
2211           + 0.5 * w * mdo.RampWearCostCoeff(:,1) .* (mdo.flow(1,j,1).mpc.gen(:,PG) - mdo.InitialPg).^2;
2212     end
2213     % Then the remaining periods
2214     for t = 2:nt
2215       for j2 = 1:mdo.idx.nj(t)
2216         for j1 = 1:mdo.idx.nj(t-1)
2217           w = mdo.tstep(t).TransMat(j2,j1) * mdo.CostWeights(1, j1, t-1);
2218           mdo.results.ExpectedRampCost(:, t) = mdo.results.ExpectedRampCost(:, t) ...
2219               + 0.5 * w * mdo.RampWearCostCoeff(:,t) .* (mdo.flow(t,j2,1).mpc.gen(:,PG) - mdo.flow(t-1,j1,1).mpc.gen(:,PG)) .^2;
2220         end
2221       end
2222     end
2223     % Finally, if there is a terminal state problem, apply cost to
2224     if ~mdo.OpenEnded
2225       for j = 1:mdi.idx.nj(nt)
2226         w = mdi.tstep(t+1).TransMat(1, j) * mdi.CostWeights(1, j, nt);
2227         mdo.results.ExpectedRampCost(:, nt+1) = 0.5 * w * mdo.RampWearCostCoeff(:,nt+1) .* (mdo.TerminalPg - mdo.flow(nt,j,1).mpc.gen(:,PG)) .^2;
2228       end
2229     end
2230     % Compute expected dispatch, conditional on making it to the
2231     % corresponding period
2232     mdo.results.ExpectedDispatch = zeros(ng, nt);
2233     for t = 1:nt
2234       pp = sum(mdo.CostWeights(1,1:mdo.idx.nj(t),t)');    % gamma(t+1)
2235       for j = 1:mdo.idx.nj(t)
2236         mdo.results.ExpectedDispatch(:,t) = mdo.results.ExpectedDispatch(:,t) + ...
2237               mdo.CostWeights(1,j,t)/pp * mdo.flow(t,j,1).mpc.gen(:,PG);
2238       end
2239     end
2240     % If Cyclic storage, pull InitialStorage value out of x
2241     if ns && mdo.Storage.ForceCyclicStorage
2242       mdo.Storage.InitialStorage = baseMVA * om.get_soln('var', 'S0');
2243     end
2244     % Compute expected storage state trajectory
2245     mdo.Storage.ExpectedStorageState = zeros(ns,nt);
2246     if ns
2247       for t = 1:nt
2248         pp = sum(mdo.CostWeights(1,1:mdo.idx.nj(t),t)');    %% gamma(t+1)
2249         for j = 1:mdo.idx.nj(t)
2250           Lfj = mdo.tstep(t).Lf( j:mdo.idx.nj(t):(ns-1)*mdo.idx.nj(t)+j, :);
2251           Ngj = mdo.tstep(t).Ng( j:mdo.idx.nj(t):(ns-1)*mdo.idx.nj(t)+j, :);
2252           Nhj = mdo.tstep(t).Nh( j:mdo.idx.nj(t):(ns-1)*mdo.idx.nj(t)+j, :);
2253           mdo.Storage.ExpectedStorageState(:,t) = ...
2254               mdo.Storage.ExpectedStorageState(:,t) + ...
2255                   baseMVA * mdo.CostWeights(1,j,t)/pp * ...
2256                   ( Lfj * mdo.Storage.InitialStorage/baseMVA + ...
2257                     (Ngj + Nhj) * mdo.QP.x );
2258         end
2259       end
2260       mdo.Storage.ExpectedStorageDispatch = ...
2261           mdo.results.ExpectedDispatch(mdo.Storage.UnitIdx, :);
2262     end
2263     % If there is a dynamical system, extract the state vectors and outputs
2264     % from the solution
2265     if ntds
2266       if nzds
2267         mdo.results.Z = zeros(nzds, ntds);
2268         for t = 1:ntds
2269           mdo.results.Z(:,t) = om.get_soln('var', 'Z', {t});
2270         end
2271       end
2272       mdo.results.Y = zeros(nyds, ntds);
2273       if nyds
2274         for t = 1:ntds
2275           AA = om.params_lin_constraint('DSy', {t});
2276           mdo.results.Y(:, t) = AA * mdo.QP.x;
2277         end
2278       end
2279     end
2280     mdo.results.f = mdo.QP.f;
2281   end   % if success
2282   mdo.results.success = success;
2283   mdo.results.SolveTime = toc(t0);
2284 end     % if mpopt.most.solve_model
2285 
2286 mdo.results.SetupTime = et_setup;
2287 
2288 if verbose
2289   fprintf('- MOST: Done.\n\n');
2290 end