most

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

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

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

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