0001 function mdo = most(mdi, mpopt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064 tmptime(1,:) = clock;
0065
0066
0067 if nargin < 2
0068 mpopt = mpoption;
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
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
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
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
0206 MinStorageLevel = ones(ns, 1) * MinStorageLevel;
0207 end
0208 if size(MinStorageLevel, 2) == 1 && nt > 1
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
0217 MaxStorageLevel = ones(ns, 1) * MaxStorageLevel;
0218 end
0219 if size(MaxStorageLevel, 2) == 1 && nt > 1
0220 MaxStorageLevel = MaxStorageLevel * ones(1, nt);
0221 end
0222 if isempty(mdi.Storage.InEff)
0223 InEff = 1;
0224 else
0225 InEff = mdi.Storage.InEff;
0226 end
0227 if size(InEff, 1) == 1 && ns > 1
0228 InEff = ones(ns, 1) * InEff;
0229 end
0230 if size(InEff, 2) == 1 && nt > 1
0231 InEff = InEff * ones(1, nt);
0232 end
0233 if isempty(mdi.Storage.OutEff)
0234 OutEff = 1;
0235 else
0236 OutEff = mdi.Storage.OutEff;
0237 end
0238 if size(OutEff, 1) == 1 && ns > 1
0239 OutEff = ones(ns, 1) * OutEff;
0240 end
0241 if size(OutEff, 2) == 1 && nt > 1
0242 OutEff = OutEff * ones(1, nt);
0243 end
0244 if isempty(mdi.Storage.LossFactor)
0245 LossFactor = 0;
0246 else
0247 LossFactor = mdi.Storage.LossFactor;
0248 end
0249 if size(LossFactor, 1) == 1 && ns > 1
0250 LossFactor = ones(ns, 1) * LossFactor;
0251 end
0252 if size(LossFactor, 2) == 1 && nt > 1
0253 LossFactor = LossFactor * ones(1, nt);
0254 end
0255 if isempty(mdi.Storage.rho)
0256 rho = 1;
0257 else
0258 rho = mdi.Storage.rho;
0259 end
0260 if size(rho, 1) == 1 && ns > 1
0261 rho = ones(ns, 1) * rho;
0262 end
0263 if size(rho, 2) == 1 && nt > 1
0264 rho = rho * ones(1, nt);
0265 end
0266 if isempty(mdi.Storage.InitialStorageLowerBound)
0267 if mo.ForceCyclicStorage
0268 mdi.Storage.InitialStorageLowerBound = MinStorageLevel(:, 1);
0269 else
0270 mdi.Storage.InitialStorageLowerBound = mdi.Storage.InitialStorage;
0271 end
0272 end
0273 if isempty(mdi.Storage.InitialStorageUpperBound)
0274 if mo.ForceCyclicStorage
0275 mdi.Storage.InitialStorageUpperBound = MaxStorageLevel(:, 1);
0276 else
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);
0308 end
0309 mdi.idx.ntds = ntds;
0310 mdi.idx.nzds = nzds;
0311 mdi.idx.nyds = nyds;
0312
0313
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
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
0337
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
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
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
0368
0369
0370
0371 mdi.UC.CommitSched = (mdi.UC.CommitKey >= 0);
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);
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);
0381 if nn > 0
0382 mdi.UC.CommitSched(i, 1:nn) = 1;
0383 end
0384 end
0385 end
0386 end
0387 end
0388
0389
0390
0391
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
0398
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
0415
0416
0417
0418 for t = 1:nt
0419
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
0424 if t == 1
0425 scenario_probs = mdi.tstep(1).TransMat;
0426 else
0427 scenario_probs = mdi.tstep(t).TransMat * mdi.CostWeights(1, 1:mdi.idx.nj(t-1), t-1)';
0428 end
0429 mdi.StepProb(t) = sum(scenario_probs);
0430 if mdi.SecurityConstrained
0431 for j = 1:mdi.idx.nj(t)
0432 [tmp, ii] = sort(mdi.cont(t,j).contab(:, CT_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 ...
0438 && mdi.flow(t,j,1).mpc.gen(contab(l, CT_ROW), GEN_STATUS) <= 0
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 ...
0442 && mdi.flow(t,j,1).mpc.branch(contab(l, CT_ROW), BR_STATUS) <= 0
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
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
0486
0487
0488
0489
0490
0491 if UC
0492 if ~isfield(mdi.UC, 'c00') || isempty(mdi.UC.c00)
0493 mdi.UC.c00 = zeros(ng, nt);
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
0518
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
0533
0534
0535
0536
0537
0538
0539
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
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
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
0583 pmin = genmask .* (min(mpc.gen(:, PMIN) / baseMVA, 0) - 1);
0584 pmax = genmask .* (max(mpc.gen(:, PMAX) / baseMVA, 0) + 1);
0585 else
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);
0596 if nrz > 1
0597 r.rgens = any(r.zones);
0598 else
0599 r.rgens = r.zones;
0600 end
0601 r.igr = find(r.rgens);
0602 ngr = length(r.igr);
0603
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
0614 if size(r.cost, 1) < ng
0615 r.original.cost = r.cost;
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;
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;
0632 Rmax = Inf(ngr, 1);
0633 kk = find(mpc.gen(r.igr, RAMP_10));
0634 Rmax(kk) = mpc.gen(r.igr(kk), RAMP_10);
0635 kk = find(r.qty(r.igr) < Rmax);
0636 Rmax(kk) = r.qty(r.igr(kk));
0637 Rmax = Rmax / baseMVA;
0638 om = add_vars(om, 'R', {t,j,k}, ngr, [], zeros(ngr, 1), Rmax);
0639 end
0640 end
0641
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
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
0650
0651
0652
0653
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
0658 end
0659
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
0666
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
0679
0680
0681
0682
0683
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
0696
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
0704
0705
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
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
0728
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
0735 end
0736 for t = 1:nt
0737 om = add_vars(om, 'Sm', {t}, ns, [], MinStorageLevel(:,t)/baseMVA, []);
0738 end
0739 end
0740
0741
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
0748
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
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));
0769 for t = 1:nt
0770 umin = zeros(ng, 1);
0771 umax = ones(ng, 1);
0772
0773 if ~mdi.UC.CyclicCommitment
0774
0775 umin( (mdi.UC.InitialState > 0) & ...
0776 (t+mdi.UC.InitialState-mdi.UC.MinUp <= 0) ) = 1;
0777
0778 umax( (mdi.UC.InitialState < 0) & ...
0779 (t-mdi.UC.InitialState-mdi.UC.MinDown <= 0) ) = 0;
0780 end
0781
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
0788 vt = vt0;
0789 vt(umin == umax) = 'C';
0790
0791 om = add_vars(om, 'u', {t}, ng, zeros(ng, 1), umin, umax, vt);
0792 end
0793
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
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
0803
0804
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
0814 qmin = genmask .* (min(mpc.gen(:, QMIN) / baseMVA, 0) - 1);
0815 qmax = genmask .* (max(mpc.gen(:, QMAX) / baseMVA, 0) + 1);
0816 else
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
0829
0830
0831
0832
0833
0834
0835
0836
0837
0838
0839
0840
0841
0842
0843
0844
0845
0846
0847
0848
0849
0850
0851
0852
0853
0854
0855
0856
0857
0858
0859
0860
0861
0862
0863
0864
0865
0866
0867 if verbose
0868 fprintf('- Building expected storage-tracking mechanism.\n');
0869 end
0870 if ns
0871
0872 vv = get_idx(om);
0873 for t = 1:nt
0874 nsxnjt = ns*mdi.idx.nj(t);
0875
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
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);
0898 mdi.tstep(t).Mh = sparse(nsxnjt, nvars);
0899 mdi.tstep(t).Ng = B2 * G;
0900 mdi.tstep(t).Nh = B2 * H;
0901 else
0902
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);
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
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
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
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
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;
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
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
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
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
1014
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
1031
1032
1033
1034
1035
1036
1037
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
1054
1055
1056
1057
1058
1059
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
1077
1078
1079
1080
1081
1082
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
1097
1098
1099
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
1114
1115
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
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
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
1148
1149 for t = 1:nt-1
1150 for j1 = 1:mdi.idx.nj(t)
1151 for j2 = 1:mdi.idx.nj(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
1163
1164
1165
1166
1167 if ~mdi.OpenEnded
1168
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
1178
1179
1180 om = add_constraints(om, 'Rrm', {nt, nj_max, nj_max});
1181
1182
1183 for t = 1:nt-1
1184 for j1 = 1:mdi.idx.nj(t)
1185 for j2 = 1:mdi.idx.nj(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
1197
1198
1199
1200
1201 if ~mdi.OpenEnded
1202
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
1214 if ns
1215 if verbose
1216 fprintf(' - Building storage constraints.\n');
1217 end
1218
1219
1220 om = add_constraints(om, 'Sm', {nt, nj_max});
1221 if mdi.Storage.ForceCyclicStorage
1222
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
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
1239
1240
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
1264 om = add_constraints(om, 'Sp', {nt, nj_max});
1265 if mdi.Storage.ForceCyclicStorage
1266
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
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
1283
1284
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
1308
1309
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
1313 if mdi.Storage.ForceCyclicStorage
1314
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
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
1330
1331
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
1357
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
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
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
1378
1379
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
1407
1408 if mdi.Storage.ForceExpectedTerminalStorage && mdi.Storage.ForceCyclicStorage
1409 error('most: ForceExpectedTerminalStorage and ForceCyclicStorage cannot be simultaneously true.');
1410 end
1411 if ns
1412
1413 if mdi.Storage.ForceExpectedTerminalStorage
1414
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
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
1452 if nzds || nyds
1453 if verbose
1454 fprintf(' - Building dynamical system constraints.\n');
1455 end
1456
1457
1458
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
1471
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
1477
1478 A = mdi.dstep(t).B * mdi.tstep(t).E;
1479 else
1480
1481
1482
1483
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
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
1511 if UC
1512 if verbose
1513 fprintf(' - Building unit commitment constraints.\n');
1514 end
1515
1516 om = add_constraints(om, 'uvw', {nt});
1517 for t = 1:nt
1518 if t == 1
1519
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
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
1538
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
1544 for tt = 1:length(ti)
1545 if ti(tt) < 1
1546 ti(tt) = nt + ti(tt);
1547 end
1548 end
1549 end
1550
1551
1552
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
1565
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
1571 for tt = 1:length(ti)
1572 if ti(tt) < 1
1573 ti(tt) = nt + ti(tt);
1574 end
1575 end
1576 end
1577
1578
1579
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
1592
1593
1594
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
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
1627
1628
1629
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
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
1669
1670
1671
1672
1673
1674
1675 c1 = 0;
1676
1677
1678
1679
1680
1681
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);
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
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
1708
1709
1710
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
1724
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);
1738
1739
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)
1745
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
1759
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
1781
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
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
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
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
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
1825 if ~mdi.OpenEnded
1826
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
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
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
1867 vv = get_idx(om);
1868 for t = 1:nt
1869
1870
1871 for j = 1:mdi.idx.nj(t)
1872
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
1908
1909
1910
1911 Cfstor(vv.i1.S0:vv.iN.S0) = ...
1912 Cfstor(vv.i1.S0:vv.iN.S0) + ...
1913 baseMVA * mdi.Storage.InitialStorageCost';
1914
1915
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
1924
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
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
1946
1947
1948
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
1966
1967
1968
1969
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
1984
1985
1986
1987
1988
1989
1990 tmptime(2,:) = clock;
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
2030
2031
2032 mdo = mdi;
2033 if mpopt.most.solve_model
2034
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
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
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;
2104
2105 if mdo.DCMODEL
2106 mpc.bus(:, VM) = 1;
2107 end
2108
2109 mpc.gen(:, PG) = baseMVA * mdo.QP.x(vv.i1.Pg(t,j,k):vv.iN.Pg(t,j,k));
2110
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
2119 mpc.bus(:, VA) = (180/pi) * mdo.QP.x(vv.i1.Va(t,j,k):vv.iN.Va(t,j,k));
2120
2121
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
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
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
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
2148
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;
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;
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;
2199 end
2200 end
2201
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
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
2216
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
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
2244 mdo.results.RrpPrices = zeros(ng, mdo.idx.ntramp);
2245 mdo.results.RrmPrices = zeros(ng, mdo.idx.ntramp);
2246
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
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
2265 mdo.results.ExpectedRampCost = zeros(ng, mdo.idx.ntramp+1);
2266
2267 for j = 1:mdi.idx.nj(1)
2268 w = mdo.tstep(1).TransMat(j,1);
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
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
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
2290
2291 mdo.results.ExpectedDispatch = zeros(ng, nt);
2292 for t = 1:nt
2293 pp = sum(mdo.CostWeights(1,1:mdo.idx.nj(t),t)');
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
2300 if ns && mdo.Storage.ForceCyclicStorage
2301 mdo.Storage.InitialStorage = baseMVA * mdo.QP.x(vv.i1.S0:vv.iN.S0);
2302 end
2303
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)');
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
2323
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
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;