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
0065
0066
0067
0068
0069 t0 = tic;
0070
0071
0072 if nargin < 2
0073 mpopt = mpoption;
0074 end
0075
0076 verbose = mpopt.verbose;
0077
0078 if verbose
0079 fprintf('\n=============================================================================\n');
0080 fprintf( ' MATPOWER Optimal Scheduling Tool -- MOST Version %s\n', mostver());
0081 fprintf( ' A multiperiod stochastic secure OPF with unit commitment\n');
0082 fprintf( ' ----- Built on MATPOWER -----\n');
0083 fprintf( ' by Carlos E. Murillo-Sanchez, Universidad Nacional de Colombia--Manizales\n');
0084 fprintf( ' and Ray D. Zimmerman, Cornell University\n');
0085 fprintf( ' (c) 2012-2020 Power Systems Engineering Research Center (PSERC) \n');
0086 fprintf( '=============================================================================\n');
0087 end
0088
0089
0090 if mpopt.most.solve_model && ~mpopt.most.resolve_new_cost
0091 mpopt = mpoption(mpopt, 'most.build_model', 1);
0092 end
0093 if ~mpopt.most.build_model && ~mpopt.most.solve_model
0094 error('most: Ah ... are you sure you want to do nothing? (either ''most.build_model'' or ''most.solve_model'' must be true)');
0095 end
0096
0097
0098 ng = size(mdi.mpc.gen, 1);
0099 nt = mdi.idx.nt;
0100 ns = length(mdi.Storage.UnitIdx);
0101 mdi.idx.ng = ng;
0102 mdi.idx.ns = ns;
0103 baseMVA = mdi.mpc.baseMVA;
0104 for t = 1:nt
0105 mdi.idx.nj(t) = length(mdi.tstep(t).OpCondSched);
0106 end
0107 if ~isfield(mdi, 'UC') || ~isfield(mdi.UC, 'CommitSched') || ...
0108 isempty(mdi.UC.CommitSched)
0109 if isfield(mdi, 'CommitSched') && ~isempty(mdi.CommitSched)
0110 warning('----- most: MDI.CommitSched has moved to MDI.UC.CommitSched, please update your code. -----');
0111 mdi.UC.CommitSched = mdi.CommitSched;
0112 else
0113 error('most: commitment schedule must be provided in MSPD_IN.UC.CommitSched');
0114 end
0115 end
0116
0117
0118 UC = mpopt.most.uc.run;
0119 if UC == -1
0120 if isempty(mdi.UC.CommitKey)
0121 UC = 0;
0122 else
0123 UC = 1;
0124 end
0125 end
0126 mo = struct(...
0127 'DCMODEL', mpopt.most.dc_model, ...
0128 'IncludeFixedReserves', mpopt.most.fixed_res, ...
0129 'SecurityConstrained', mpopt.most.security_constraints, ...
0130 'QCoordination', mpopt.most.q_coordination, ...
0131 'ForceCyclicStorage', mpopt.most.storage.cyclic, ...
0132 'CyclicCommitment', mpopt.most.uc.cyclic, ...
0133 'ForceExpectedTerminalStorage', mpopt.most.storage.terminal_target, ...
0134 'alpha', mpopt.most.alpha ...
0135 );
0136 if mo.IncludeFixedReserves == -1
0137 if isfield(mdi, 'FixedReserves') && isfield(mdi.FixedReserves(1,1,1), 'req')
0138 mo.IncludeFixedReserves = 1;
0139 else
0140 mo.IncludeFixedReserves = 0;
0141 end
0142 end
0143 if mo.SecurityConstrained == -1
0144 if isfield(mdi, 'cont') && isfield(mdi.cont(1,1), 'contab') && ...
0145 ~isempty(mdi.cont(1,1).contab)
0146 mo.SecurityConstrained = 1;
0147 else
0148 mo.SecurityConstrained = 0;
0149 end
0150 end
0151 if mo.ForceExpectedTerminalStorage == -1
0152 if isempty(mdi.Storage.ExpectedTerminalStorageAim) && ...
0153 isempty(mdi.Storage.ExpectedTerminalStorageMin) && ...
0154 isempty(mdi.Storage.ExpectedTerminalStorageMax)
0155 mo.ForceExpectedTerminalStorage = 0;
0156 else
0157 mo.ForceExpectedTerminalStorage = 1;
0158 end
0159 end
0160 if mo.IncludeFixedReserves && ~(isfield(mdi, 'FixedReserves') && ...
0161 isfield(mdi.FixedReserves(1,1,1), 'req'))
0162 error('most: MDI.FixedReserves(t,j,k) must be specified when MPOPT.most.fixed_res = 1');
0163 end
0164 if mo.SecurityConstrained && ~(isfield(mdi, 'cont') && ...
0165 isfield(mdi.cont(1,1), 'contab') && ~isempty(mdi.cont(1,1).contab))
0166 error('most: MDI.cont(t,j).contab cannot be empty when MPOPT.most.security_constraints = 1');
0167 end
0168 if mo.IncludeFixedReserves && mo.SecurityConstrained
0169 warning('most: Using MPOPT.most.fixed_res = 1 and MPOPT.most.security_constraints = 1 together is not recommended.');
0170 end
0171 if mo.ForceExpectedTerminalStorage == 1;
0172 if mo.ForceCyclicStorage
0173 error('most: storage model cannot be both cyclic and include a terminal target value; must change MPOPT.most.storage.cyclic or MPOPT.most.storage.terminal_target');
0174 end
0175 if ns && isempty(mdi.Storage.ExpectedTerminalStorageAim) && ...
0176 isempty(mdi.Storage.ExpectedTerminalStorageMin) && ...
0177 isempty(mdi.Storage.ExpectedTerminalStorageMax)
0178 error('most: MDI.Storage.ExpectedTerminalStorageAim|Min|Max cannot all be empty when MPOPT.most.storage.terminal_target = 1');
0179 end
0180 if ~isempty(mdi.Storage.ExpectedTerminalStorageAim)
0181 mdi.Storage.ExpectedTerminalStorageMin = mdi.Storage.ExpectedTerminalStorageAim;
0182 mdi.Storage.ExpectedTerminalStorageMax = mdi.Storage.ExpectedTerminalStorageAim;
0183 end
0184 end
0185 if UC && (~isfield(mdi.UC, 'CommitKey') || isempty(mdi.UC.CommitKey))
0186 error('most: cannot run unit commitment without specifying MDI.UC.CommitKey');
0187 end
0188
0189 if ns
0190 if isempty(mdi.Storage.InitialStorage)
0191 error('most: Storage.InitialStorage must be specified');
0192 end
0193 if isempty(mdi.Storage.TerminalChargingPrice0)
0194 mdi.Storage.TerminalChargingPrice0 = mdi.Storage.TerminalStoragePrice;
0195 end
0196 if isempty(mdi.Storage.TerminalDischargingPrice0)
0197 mdi.Storage.TerminalDischargingPrice0 = mdi.Storage.TerminalStoragePrice;
0198 end
0199 if isempty(mdi.Storage.TerminalChargingPriceK)
0200 mdi.Storage.TerminalChargingPriceK = mdi.Storage.TerminalStoragePrice;
0201 end
0202 if isempty(mdi.Storage.TerminalDischargingPriceK)
0203 mdi.Storage.TerminalDischargingPriceK = mdi.Storage.TerminalStoragePrice;
0204 end
0205 if isempty(mdi.Storage.MinStorageLevel)
0206 error('most: Storage.MinStorageLevel must be specified');
0207 else
0208 MinStorageLevel = mdi.Storage.MinStorageLevel;
0209 end
0210 if size(MinStorageLevel, 1) == 1 && ns > 1
0211 MinStorageLevel = ones(ns, 1) * MinStorageLevel;
0212 end
0213 if size(MinStorageLevel, 2) == 1 && nt > 1
0214 MinStorageLevel = MinStorageLevel * ones(1, nt);
0215 end
0216 if isempty(mdi.Storage.MaxStorageLevel)
0217 error('most: Storage.MaxStorageLevel must be specified');
0218 else
0219 MaxStorageLevel = mdi.Storage.MaxStorageLevel;
0220 end
0221 if size(MaxStorageLevel, 1) == 1 && ns > 1
0222 MaxStorageLevel = ones(ns, 1) * MaxStorageLevel;
0223 end
0224 if size(MaxStorageLevel, 2) == 1 && nt > 1
0225 MaxStorageLevel = MaxStorageLevel * ones(1, nt);
0226 end
0227 if isempty(mdi.Storage.InEff)
0228 InEff = 1;
0229 else
0230 InEff = mdi.Storage.InEff;
0231 end
0232 if size(InEff, 1) == 1 && ns > 1
0233 InEff = ones(ns, 1) * InEff;
0234 end
0235 if size(InEff, 2) == 1 && nt > 1
0236 InEff = InEff * ones(1, nt);
0237 end
0238 if isempty(mdi.Storage.OutEff)
0239 OutEff = 1;
0240 else
0241 OutEff = mdi.Storage.OutEff;
0242 end
0243 if size(OutEff, 1) == 1 && ns > 1
0244 OutEff = ones(ns, 1) * OutEff;
0245 end
0246 if size(OutEff, 2) == 1 && nt > 1
0247 OutEff = OutEff * ones(1, nt);
0248 end
0249 if isempty(mdi.Storage.LossFactor)
0250 LossFactor = 0;
0251 else
0252 LossFactor = mdi.Storage.LossFactor;
0253 end
0254 if size(LossFactor, 1) == 1 && ns > 1
0255 LossFactor = ones(ns, 1) * LossFactor;
0256 end
0257 if size(LossFactor, 2) == 1 && nt > 1
0258 LossFactor = LossFactor * ones(1, nt);
0259 end
0260 if isempty(mdi.Storage.rho)
0261 rho = 1;
0262 else
0263 rho = mdi.Storage.rho;
0264 end
0265 if size(rho, 1) == 1 && ns > 1
0266 rho = ones(ns, 1) * rho;
0267 end
0268 if size(rho, 2) == 1 && nt > 1
0269 rho = rho * ones(1, nt);
0270 end
0271 if isempty(mdi.Storage.InitialStorageLowerBound)
0272 if mo.ForceCyclicStorage
0273 mdi.Storage.InitialStorageLowerBound = MinStorageLevel(:, 1);
0274 else
0275 mdi.Storage.InitialStorageLowerBound = mdi.Storage.InitialStorage;
0276 end
0277 end
0278 if isempty(mdi.Storage.InitialStorageUpperBound)
0279 if mo.ForceCyclicStorage
0280 mdi.Storage.InitialStorageUpperBound = MaxStorageLevel(:, 1);
0281 else
0282 mdi.Storage.InitialStorageUpperBound = mdi.Storage.InitialStorage;
0283 end
0284 end
0285
0286 LossCoeff = mdi.Delta_T * LossFactor/2;
0287 beta1 = (1-LossCoeff) ./ (1+LossCoeff);
0288 beta2 = 1 ./ (1+LossCoeff);
0289 beta3 = 1 ./ (1/(1-mo.alpha) + LossCoeff);
0290 beta4 = mo.alpha/(1-mo.alpha) * beta2 .* beta3;
0291 beta5 = beta1 ./ beta2 .* (beta3 + beta4);
0292 beta2EtaIn = mdi.Delta_T * beta2 .* InEff;
0293 beta2overEtaOut = mdi.Delta_T * beta2 ./ OutEff;
0294 beta3EtaIn = mdi.Delta_T * beta3 .* InEff;
0295 beta3overEtaOut = mdi.Delta_T * beta3 ./ OutEff;
0296 beta4EtaIn = mdi.Delta_T * beta4 .* InEff;
0297 beta4overEtaOut = mdi.Delta_T * beta4 ./ OutEff;
0298 diagBeta2EtaIn1 = spdiags(beta2EtaIn(:,1), 0, ns, ns);
0299 diagBeta2overEtaOut1 = spdiags(beta2overEtaOut(:,1), 0, ns, ns);
0300 diagBeta3EtaIn1 = spdiags(beta3EtaIn(:,1), 0, ns, ns);
0301 diagBeta3overEtaOut1 = spdiags(beta3overEtaOut(:,1), 0, ns, ns);
0302 diagBeta4EtaIn1 = spdiags(beta4EtaIn(:,1), 0, ns, ns);
0303 diagBeta4overEtaOut1 = spdiags(beta4overEtaOut(:,1), 0, ns, ns);
0304 end
0305 if ~isfield(mdi.idx, 'ntds') || isempty(mdi.idx.ntds) || ~mdi.idx.ntds
0306 ntds = 0;
0307 nzds = 0;
0308 nyds = 0;
0309 else
0310 ntds = mdi.idx.ntds;
0311 nzds = size(mdi.dstep(1).A, 1);
0312 nyds = size(mdi.dstep(1).D, 1);
0313 end
0314 mdi.idx.ntds = ntds;
0315 mdi.idx.nzds = nzds;
0316 mdi.idx.nyds = nyds;
0317
0318
0319 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0320 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0321 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0322 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0323 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0324 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0325 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0326 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0327 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0328 [CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
0329 CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
0330 CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
0331 CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
0332 CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
0333 CT_MODCOST_X] = idx_ct;
0334
0335
0336 nb = size(mdi.mpc.bus, 1);
0337 if any(mdi.mpc.bus(:, BUS_I) ~= (1:nb)')
0338 error('most: buses must be numbered consecutively in bus matrix; use ext2int() to convert to internal ordering')
0339 end
0340
0341
0342
0343 mdi.mpc.bus(:, MU_VMIN) = 0;
0344 mdi.mpc.gen(:, MU_QMIN) = 0;
0345 mdi.mpc.branch(:,MU_ANGMAX) = 0;
0346 mdi.mpc.f = 0;
0347 mdi.mpc.et = 0;
0348 mdi.mpc.success = 1;
0349
0350 if mpopt.most.build_model
0351 if verbose
0352 fprintf('- Building indexing structures.\n');
0353 end
0354
0355
0356 mdi.DCMODEL = mo.DCMODEL;
0357 mdi.IncludeFixedReserves = mo.IncludeFixedReserves;
0358 mdi.SecurityConstrained = mo.SecurityConstrained;
0359 mdi.QCoordination = mo.QCoordination;
0360 mdi.Storage.ForceCyclicStorage = mo.ForceCyclicStorage;
0361 mdi.Storage.ForceExpectedTerminalStorage = mo.ForceExpectedTerminalStorage;
0362 mdi.UC.run = UC;
0363 mdi.UC.CyclicCommitment = mo.CyclicCommitment;
0364 mdi.alpha = mo.alpha;
0365 if ~isfield(mdi, 'OpenEnded'), mdi.OpenEnded = 1; end
0366
0367 if UC
0368
0369 if any(mdi.UC.MinUp < 1) && any(mdi.UC.MinUp < 1)
0370 error('most: UC.MinUp and UC.MinDown must all be >= 1');
0371 end
0372
0373
0374
0375
0376 mdi.UC.CommitSched = (mdi.UC.CommitKey >= 0);
0377 if ~mdi.UC.CyclicCommitment
0378 for i = 1:ng
0379 if mdi.UC.InitialState(i) < 0
0380 nn = mdi.UC.MinDown(i) + mdi.UC.InitialState(i);
0381 if nn > 0
0382 mdi.UC.CommitSched(i, 1:nn) = 0;
0383 end
0384 elseif mdi.UC.InitialState(i) > 0
0385 nn = mdi.UC.MinUp(i) - mdi.UC.InitialState(i);
0386 if nn > 0
0387 mdi.UC.CommitSched(i, 1:nn) = 1;
0388 end
0389 end
0390 end
0391 end
0392 end
0393
0394
0395
0396
0397 for t = 1:nt
0398 for j = 1:mdi.idx.nj(t)
0399 mpc = mdi.mpc;
0400 mpc.gen(:, GEN_STATUS) = mdi.UC.CommitSched(:, t);
0401
0402
0403
0404 if isfield(mdi.offer(t), 'gencost') && ~isempty(mdi.offer(t).gencost)
0405 mpc.gencost = mdi.offer(t).gencost;
0406 end
0407
0408 if ~isempty(mdi.tstep(t).OpCondSched(j).tab)
0409 changelist = unique(mdi.tstep(t).OpCondSched(j).tab(:, CT_LABEL));
0410 for label = changelist'
0411 mpc = apply_changes(label, mpc, mdi.tstep(t).OpCondSched(j).tab);
0412 end
0413 end
0414 mdi.flow(t,j,1).mpc = mpc;
0415 mdi.idx.nb(t,j,1) = size(mdi.flow(t,j,1).mpc.bus, 1);
0416 mdi.idx.ny(t,j,1) = length(find(mdi.flow(t,j,1).mpc.gencost(:, MODEL) == PW_LINEAR));
0417 end
0418 end
0419
0420
0421
0422
0423 for t = 1:nt
0424
0425 if ~isfield(mdi.tstep(t), 'TransMask') || isempty(mdi.tstep(t).TransMask)
0426 mdi.tstep(t).TransMask = ones(size(mdi.tstep(t).TransMat));
0427 end
0428
0429 if t == 1
0430 scenario_probs = mdi.tstep(1).TransMat;
0431 else
0432 scenario_probs = mdi.tstep(t).TransMat * mdi.CostWeights(1, 1:mdi.idx.nj(t-1), t-1)';
0433 end
0434 mdi.StepProb(t) = sum(scenario_probs);
0435 if mdi.SecurityConstrained
0436 for j = 1:mdi.idx.nj(t)
0437 [tmp, ii] = sort(mdi.cont(t,j).contab(:, CT_LABEL));
0438 contab = mdi.cont(t,j).contab(ii, :);
0439 rowdecomlist = ones(size(contab,1), 1);
0440 for l = 1:size(contab, 1)
0441 if contab(l, CT_TABLE) == CT_TGEN && contab(l, CT_COL) == GEN_STATUS ...
0442 && contab(l, CT_CHGTYPE) == CT_REP && contab(l, CT_NEWVAL) == 0 ...
0443 && mdi.flow(t,j,1).mpc.gen(contab(l, CT_ROW), GEN_STATUS) <= 0
0444 rowdecomlist(l) = 0;
0445 elseif contab(l, CT_TABLE) == CT_TBRCH && contab(l, CT_COL) == BR_STATUS ...
0446 && contab(l, CT_CHGTYPE) == CT_REP && contab(l, CT_NEWVAL) == 0 ...
0447 && mdi.flow(t,j,1).mpc.branch(contab(l, CT_ROW), BR_STATUS) <= 0
0448 rowdecomlist(l) = 0;
0449 end
0450 end
0451 contab = contab(rowdecomlist ~= 0, :);
0452 mdi.cont(t, j).contab = contab;
0453 clist = unique(contab(:, CT_LABEL));
0454 mdi.idx.nc(t, j) = length(clist);
0455 k = 2;
0456 for label = clist'
0457 mdi.flow(t, j, k).mpc = apply_changes(label, mdi.flow(t, j, 1).mpc, contab);
0458 ii = find( label == contab(:, CT_LABEL) );
0459 mdi.CostWeights(k, j, t) = contab(ii(1), CT_PROB);
0460 mdi.idx.nb(t, j, k) = size(mdi.flow(t, j, k).mpc.bus, 1);
0461 mdi.idx.ny(t, j, k) = length(find(mdi.flow(t, j, 1).mpc.gencost(:, MODEL) == PW_LINEAR));
0462 k = k + 1;
0463 end
0464 mdi.CostWeights(1, j, t) = 1 - sum(mdi.CostWeights(2:mdi.idx.nc(t,j)+1, j, t));
0465 mdi.CostWeights(1:mdi.idx.nc(t,j)+1, j, t) = scenario_probs(j) * mdi.CostWeights(1:mdi.idx.nc(t,j)+1, j, t);
0466 end
0467 else
0468 for j = 1:mdi.idx.nj(t)
0469 mdi.idx.nc(t, j) = 0;
0470 mdi.CostWeights(1, j, t) = scenario_probs(j);
0471 end
0472 end
0473 end
0474
0475
0476 if mdi.SecurityConstrained && mdi.alpha ~= 0
0477 for t = 1:nt
0478 for j = 1:mdi.idx.nj(t)
0479 mdi.CostWeightsAdj(1, j, t) = mdi.CostWeights(1, j, t);
0480 for k = 2:mdi.idx.nc(t,j)+1
0481 mdi.CostWeightsAdj(k, j, t) = (1-mdi.alpha) * mdi.CostWeights(k, j, t);
0482 mdi.CostWeightsAdj(1, j, t) = mdi.CostWeightsAdj(1, j, t) + mdi.alpha * mdi.CostWeights(k, j, t);
0483 end
0484 end
0485 end
0486 else
0487 mdi.CostWeightsAdj = mdi.CostWeights;
0488 end
0489
0490
0491
0492
0493
0494
0495
0496 if UC
0497 if ~isfield(mdi.UC, 'c00') || isempty(mdi.UC.c00)
0498 mdi.UC.c00 = zeros(ng, nt);
0499 for t = 1:nt
0500 for j = 1:mdi.idx.nj(t)
0501 for k = 1:mdi.idx.nc(t,j)+1
0502 mpc = mdi.flow(t,j,k).mpc;
0503 c00tjk = totcost(mpc.gencost, zeros(ng,1));
0504 mdi.UC.c00(:, t) = mdi.UC.c00(:, t) + mdi.CostWeightsAdj(k, j, t) * c00tjk;
0505 c0col = COST + mpc.gencost(:,NCOST) - 1;
0506 ipoly = find(mpc.gencost(:, MODEL) == POLYNOMIAL);
0507 ipwl = find(mpc.gencost(:, MODEL) == PW_LINEAR);
0508 ii = sub2ind(size(mpc.gencost), ipoly, c0col(ipoly));
0509 mpc.gencost(ii) = mpc.gencost(ii) - c00tjk(ipoly);
0510 for i = ipwl'
0511 jj = COST+1:2:COST+2*mpc.gencost(i,NCOST)-1;
0512 mpc.gencost(i, jj) = mpc.gencost(i, jj) - c00tjk(i);
0513 end
0514 mpc.fixed_gencost = c00tjk;
0515 mdi.flow(t,j,k).mpc = mpc;
0516 end
0517 end
0518 end
0519 end
0520 end
0521
0522
0523
0524 mdi.idx.nf_total = 0;
0525 mdi.idx.nb_total = 0;
0526 for t = 1:nt
0527 for j = 1:mdi.idx.nj(t)
0528 mdi.idx.nf_total = mdi.idx.nf_total + (1 + mdi.idx.nc(t,j));
0529 for k = 1:mdi.idx.nc(t,j)+1
0530 mdi.idx.nb_total = mdi.idx.nb_total + size(mdi.flow(t, j, k).mpc.bus, 1);
0531 ii = find(mdi.flow(t,j,k).mpc.gencost(:, MODEL) == PW_LINEAR);
0532 mdi.idx.ny(t,j,k) = length(ii);
0533 end
0534 end
0535 end
0536 mdi.idx.ns_total = ns * mdi.idx.nf_total;
0537
0538
0539
0540
0541
0542
0543
0544
0545 om = opt_model;
0546 nj_max = max(mdi.idx.nj);
0547 nc_max = max(max(mdi.idx.nc));
0548 Ing = speye(ng);
0549 Ins = speye(ns);
0550 if mdi.DCMODEL
0551 om.init_indexed_name('var', 'Va', {nt, nj_max, nc_max+1});
0552 end
0553 om.init_indexed_name('var', 'Pg', {nt, nj_max, nc_max+1});
0554 om.init_indexed_name('var', 'dPp', {nt, nj_max, nc_max+1});
0555 om.init_indexed_name('var', 'dPm', {nt, nj_max, nc_max+1});
0556 om.init_indexed_name('var', 'y', {nt, nj_max, nc_max+1});
0557 if mdi.IncludeFixedReserves
0558 om.init_indexed_name('var', 'R', {nt, nj_max, nc_max+1});
0559 end
0560 for t = 1:nt
0561 for j = 1:mdi.idx.nj(t)
0562
0563 for k = 1:mdi.idx.nc(t,j)+1
0564 if mdi.DCMODEL
0565 iref = find(mdi.flow(t,j,k).mpc.bus(:,BUS_TYPE) == REF);
0566 if verbose && length(iref) > 1
0567 errstr = ['\nmost: Warning: Multiple reference buses.\n', ...
0568 ' For a system with islands, a reference bus in each island\n', ...
0569 ' may help convergence, but in a fully connected system such\n', ...
0570 ' a situation is probably not reasonable.\n\n' ];
0571 fprintf(errstr);
0572 end
0573 Va0 = mdi.flow(t,j,k).mpc.bus(:,VA)*pi/180;
0574 Va_max = Inf(mdi.idx.nb(t,j,k), 1);
0575 Va_min = -Va_max;
0576 Va_min(iref) = mdi.flow(t,j,k).mpc.bus(iref,VA)*pi/180;
0577 Va_max(iref) = Va_min(iref);
0578
0579 om.add_var('Va', {t,j,k}, mdi.idx.nb(t,j,k), Va0, Va_min, Va_max);
0580 end
0581 end
0582
0583 for k = 1:mdi.idx.nc(t,j)+1
0584 mpc = mdi.flow(t,j,k).mpc;
0585 genmask = mpc.gen(:,GEN_STATUS) > 0;
0586 p0 = genmask .* mpc.gen(:,PG) / baseMVA;
0587 if UC
0588 pmin = genmask .* (min(mpc.gen(:, PMIN) / baseMVA, 0) - 1);
0589 pmax = genmask .* (max(mpc.gen(:, PMAX) / baseMVA, 0) + 1);
0590 else
0591 pmin = genmask .* mpc.gen(:, PMIN) / baseMVA;
0592 pmax = genmask .* mpc.gen(:, PMAX) / baseMVA;
0593 end
0594 om.add_var('Pg', {t,j,k}, ng, p0, pmin, pmax);
0595 end
0596 if mdi.IncludeFixedReserves
0597 for k = 1:mdi.idx.nc(t,j)+1
0598 mpc = mdi.flow(t,j,k).mpc;
0599 r = mdi.FixedReserves(t,j,k);
0600 nrz = size(r.req, 1);
0601 if nrz > 1
0602 r.rgens = any(r.zones);
0603 else
0604 r.rgens = r.zones;
0605 end
0606 r.igr = find(r.rgens);
0607 ngr = length(r.igr);
0608
0609 if size(r.zones, 1) ~= nrz
0610 error('most: the number of rows in FixedReserves(%d,%d,%d).req (%d) and FixedReserves(%d,%d,%d).zones (%d) must match', t, j, k, nrz, t, j, k, size(r.zones, 1));
0611 end
0612 if size(r.cost, 1) ~= ng && size(r.cost, 1) ~= ngr
0613 error('most: the number of rows in FixedReserves(%d,%d,%d).cost (%d) must equal the total number of generators (%d) or the number of generators able to provide reserves (%d)', t, j, k, size(r.cost, 1), ng, ngr);
0614 end
0615 if isfield(r, 'qty') && size(r.qty, 1) ~= size(r.cost, 1)
0616 error('most: FixedReserves(%d,%d,%d).cost (%d x 1) and FixedReserves(%d,%d,%d).qty (%d x 1) must be the same dimension', t, j, k, size(r.cost, 1), t, j, k, size(r.qty, 1));
0617 end
0618
0619 if size(r.cost, 1) < ng
0620 r.original.cost = r.cost;
0621 cost = zeros(ng, 1);
0622 cost(r.igr) = r.cost;
0623 r.cost = cost;
0624 if isfield(r, 'qty')
0625 r.original.qty = r.qty;
0626 qty = zeros(ng, 1);
0627 qty(r.igr) = r.qty;
0628 r.qty = qty;
0629 end
0630 end
0631 mdi.FixedReserves(t,j,k).rgens = r.rgens;
0632 mdi.FixedReserves(t,j,k).igr = r.igr;
0633 if isfield(r, 'original')
0634 mdi.FixedReserves(t,j,k).original = r.original;
0635 end
0636 mdi.FixedReserves(t,j,k) = r;
0637 Rmax = Inf(ngr, 1);
0638 kk = find(mpc.gen(r.igr, RAMP_10));
0639 Rmax(kk) = mpc.gen(r.igr(kk), RAMP_10);
0640 kk = find(r.qty(r.igr) < Rmax);
0641 Rmax(kk) = r.qty(r.igr(kk));
0642 Rmax = Rmax / baseMVA;
0643 om.add_var('R', {t,j,k}, ngr, [], zeros(ngr, 1), Rmax);
0644 end
0645 end
0646
0647 for k = 1:mdi.idx.nc(t,j)+1
0648 om.add_var('dPp', {t,j,k}, ng, [], zeros(ng,1), []);
0649 end
0650
0651 for k = 1:mdi.idx.nc(t,j)+1
0652 om.add_var('dPm', {t,j,k}, ng, [], zeros(ng,1), []);
0653 end
0654
0655
0656
0657
0658
0659 for k = 1:mdi.idx.nc(t,j)+1
0660 om.add_var('y', {t,j,k}, mdi.idx.ny(t,j,k), [], [], []);
0661 end
0662 end
0663 end
0664
0665 om.init_indexed_name('var', 'Pc', {nt});
0666 om.init_indexed_name('var', 'Rpp', {nt});
0667 om.init_indexed_name('var', 'Rpm', {nt});
0668 for t = 1:nt
0669 om.add_var('Pc', {t}, ng);
0670
0671
0672 Rpmin = -Inf(ng,1);
0673 off = ones(ng,1);
0674 for j = 1:mdi.idx.nj(t);
0675 for k = 1:mdi.idx.nc(t,j)+1
0676 off = off & mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) <= 0;
0677 end
0678 end
0679 Rpmin(off == 1) = 0;
0680 om.add_var('Rpp', {t}, ng, [], Rpmin, mdi.offer(t).PositiveActiveReserveQuantity/baseMVA);
0681 om.add_var('Rpm', {t}, ng, [], Rpmin, mdi.offer(t).NegativeActiveReserveQuantity/baseMVA);
0682 end
0683
0684
0685
0686
0687
0688
0689 if ~mdi.OpenEnded
0690 mdi.idx.ntramp = nt;
0691 else
0692 mdi.idx.ntramp = nt - 1;
0693 end
0694 om.init_indexed_name('var', 'Rrp', {mdi.idx.ntramp});
0695 om.init_indexed_name('var', 'Rrm', {mdi.idx.ntramp});
0696 for t = 1:mdi.idx.ntramp
0697 ramp30 = mdi.flow(t,1,1).mpc.gen(:,RAMP_30)*2*mdi.Delta_T;
0698 om.add_var('Rrp', {t}, ng, [], zeros(ng,1), ...
0699 min(mdi.offer(t).PositiveLoadFollowReserveQuantity, ramp30)/baseMVA);
0700 end
0701 for t = 1:mdi.idx.ntramp
0702 ramp30 = mdi.flow(t,1,1).mpc.gen(:,RAMP_30)*2*mdi.Delta_T;
0703 om.add_var('Rrm', {t}, ng, [], zeros(ng,1), ...
0704 min(mdi.offer(t).NegativeLoadFollowReserveQuantity, ramp30)/baseMVA);
0705 end
0706
0707
0708
0709 om.init_indexed_name('var', 'Psc', {nt, nj_max, nc_max+1});
0710 om.init_indexed_name('var', 'Psd', {nt, nj_max, nc_max+1});
0711 if ns
0712 for t = 1:nt
0713 for j = 1:mdi.idx.nj(t)
0714 for k = 1:mdi.idx.nc(t,j)+1
0715 om.add_var('Psc', {t,j,k}, ns, [], [], zeros(ns,1));
0716 end
0717 end
0718 end
0719 for t = 1:nt
0720 for j = 1:mdi.idx.nj(t)
0721 for k = 1:mdi.idx.nc(t,j)+1
0722 om.add_var('Psd', {t,j,k}, ns, [], zeros(ns,1), []);
0723 end
0724 end
0725 end
0726 end
0727
0728
0729 om.init_indexed_name('var', 'Sp', {nt});
0730 om.init_indexed_name('var', 'Sm', {nt});
0731 if ns
0732 for t = 1:nt
0733 om.add_var('Sp', {t}, ns, [], [], MaxStorageLevel(:,t)/baseMVA);
0734 end
0735 for t = 1:nt
0736 om.add_var('Sm', {t}, ns, [], MinStorageLevel(:,t)/baseMVA, []);
0737 end
0738 end
0739
0740
0741 if ns && mdi.Storage.ForceCyclicStorage
0742 om.add_var('S0', ns, [], ...
0743 mdi.Storage.InitialStorageLowerBound / baseMVA, ...
0744 mdi.Storage.InitialStorageUpperBound / baseMVA);
0745 end
0746
0747
0748 if nzds
0749 om.init_indexed_name('var', 'Z', {ntds});
0750 for t = 1:ntds
0751 if t == 1
0752 zmin = mdi.z1;
0753 zmax = mdi.z1;
0754 else
0755 zmin = mdi.dstep(t).zmin;
0756 zmax = mdi.dstep(t).zmax;
0757 end
0758 z0 = (zmax - zmin) / 2;
0759 om.add_var('Z', {t}, nzds, z0, zmin, zmax);
0760 end
0761 end
0762
0763 if UC
0764 om.init_indexed_name('var', 'u', {nt});
0765 om.init_indexed_name('var', 'v', {nt});
0766 om.init_indexed_name('var', 'w', {nt});
0767 vt0 = char('B' * ones(1, ng));
0768 for t = 1:nt
0769 umin = zeros(ng, 1);
0770 umax = ones(ng, 1);
0771
0772 if ~mdi.UC.CyclicCommitment
0773
0774 umin( (mdi.UC.InitialState > 0) & ...
0775 (t+mdi.UC.InitialState-mdi.UC.MinUp <= 0) ) = 1;
0776
0777 umax( (mdi.UC.InitialState < 0) & ...
0778 (t-mdi.UC.InitialState-mdi.UC.MinDown <= 0) ) = 0;
0779 end
0780
0781 iON = find(mdi.UC.CommitKey(:,t) == 2);
0782 iOFF = find(mdi.UC.CommitKey(:,t) == -1);
0783 umin(iON) = 1;
0784 umax(iOFF) = 0;
0785
0786
0787 vt = vt0;
0788 vt(umin == umax) = 'C';
0789
0790 om.add_var('u', {t}, ng, zeros(ng, 1), umin, umax, vt);
0791 end
0792
0793 for t = 1:nt
0794 om.add_var('v', {t}, ng, zeros(ng, 1), zeros(ng, 1), ones(ng, 1));
0795 end
0796
0797 for t = 1:nt
0798 om.add_var('w', {t}, ng, zeros(ng, 1), zeros(ng, 1), ones(ng, 1));
0799 end
0800 end
0801
0802
0803
0804 if mdi.QCoordination
0805 om.init_indexed_name('var', 'Qg', {nt, nj_max, nc_max+1});
0806 for t = 1:nt
0807 for j = 1:mdi.idx.nj(t)
0808 for k = 1:mdi.idx.nc(t,j)+1
0809 mpc = mdi.flow(t,j,k).mpc;
0810 genmask = mpc.gen(:,GEN_STATUS) > 0;
0811 q0 = genmask .* mpc.gen(:, QG) / baseMVA;
0812 if UC
0813 qmin = genmask .* (min(mpc.gen(:, QMIN) / baseMVA, 0) - 1);
0814 qmax = genmask .* (max(mpc.gen(:, QMAX) / baseMVA, 0) + 1);
0815 else
0816 qmin = genmask .* mpc.gen(:, QMIN) / baseMVA;
0817 qmax = genmask .* mpc.gen(:, QMAX) / baseMVA;
0818 end
0819 om.add_var('Qg', {t,j,k}, ng, q0, qmin, qmax);
0820 end
0821 end
0822 end
0823 end
0824
0825
0826
0827 nvars = om.getN('var');
0828 mdi.idx.nvars = nvars;
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
0868
0869 if verbose
0870 fprintf('- Building expected storage-tracking mechanism.\n');
0871 end
0872 if ns
0873
0874 vv = om.get_idx();
0875 for t = 1:nt
0876 nsxnjt = ns*mdi.idx.nj(t);
0877
0878 G = sparse(nsxnjt, nvars);
0879 H = sparse(nsxnjt, nvars);
0880 B1 = sparse(nsxnjt, nsxnjt);
0881 B2 = sparse(nsxnjt, nsxnjt);
0882 for j = 1:mdi.idx.nj(t)
0883 ii = ((1:ns)'-1)*mdi.idx.nj(t)+j;
0884 jj1 = (vv.i1.Psc(t,j,1):vv.iN.Psc(t,j,1))';
0885 jj2 = (vv.i1.Psd(t,j,1):vv.iN.Psd(t,j,1))';
0886 G = G + sparse(ii, jj1, -mdi.Delta_T * InEff(:,t), nsxnjt, nvars);
0887 H = H + sparse(ii, jj2, -mdi.Delta_T ./ OutEff(:,t), nsxnjt, nvars);
0888 B1 = B1 + sparse(ii, ii, beta1(:,t), nsxnjt, nsxnjt);
0889 B2 = B2 + sparse(ii, ii, beta2(:,t), nsxnjt, nsxnjt);
0890 end
0891 if t == 1
0892
0893 jlist = [];
0894 for i=1:ns
0895 jlist = [ jlist; i*ones(mdi.idx.nj(t),1) ];
0896 end
0897 mdi.tstep(t).Li = sparse((1:nsxnjt)', jlist, 1, nsxnjt, ns);
0898 mdi.tstep(t).Lf = B1 * mdi.tstep(t).Li;
0899 mdi.tstep(t).Mg = sparse(nsxnjt, nvars);
0900 mdi.tstep(t).Mh = sparse(nsxnjt, nvars);
0901 mdi.tstep(t).Ng = B2 * G;
0902 mdi.tstep(t).Nh = B2 * H;
0903 else
0904
0905 D = sparse(nsxnjt, ns*mdi.idx.nj(t-1));
0906 p1 = mdi.CostWeights(1,1:mdi.idx.nj(t-1),t-1)';
0907 p1 = p1 / sum(p1);
0908 p2 = mdi.tstep(t).TransMat * p1;
0909 Di = spdiags(1./p2, 0, mdi.idx.nj(t), mdi.idx.nj(t)) * ...
0910 sparse(mdi.tstep(t).TransMat) * ...
0911 spdiags(p1, 0, mdi.idx.nj(t-1), mdi.idx.nj(t-1));
0912 for i = 1:ns
0913 D((i-1)*mdi.idx.nj(t)+1:i*mdi.idx.nj(t), (i-1)*mdi.idx.nj(t-1)+1:i*mdi.idx.nj(t-1)) = Di;
0914 end
0915
0916 mdi.tstep(t).Li = D * mdi.tstep(t-1).Lf;
0917 mdi.tstep(t).Lf = B1 * mdi.tstep(t).Li;
0918 mdi.tstep(t).Mg = D * mdi.tstep(t-1).Ng;
0919 mdi.tstep(t).Mh = D * mdi.tstep(t-1).Nh;
0920 mdi.tstep(t).Ng = B1 * mdi.tstep(t).Mg + B2 * G;
0921 mdi.tstep(t).Nh = B1 * mdi.tstep(t).Mh + B2 * H;
0922 end
0923 mdi.tstep(t).G = G;
0924 mdi.tstep(t).H = H;
0925 end
0926 end
0927
0928
0929 if verbose
0930 fprintf('- Building constraint submatrices.\n');
0931 end
0932 baseMVA = mdi.mpc.baseMVA;
0933 om.init_indexed_name('lin', 'Pmis', {nt, nj_max, nc_max+1});
0934 if mdi.DCMODEL
0935
0936 if verbose
0937 fprintf(' - Building DC flow constraints.\n');
0938 end
0939 om.init_indexed_name('lin', 'Pf', {nt, nj_max, nc_max+1});
0940 for t = 1:nt
0941 for j = 1:mdi.idx.nj(t)
0942 for k = 1:mdi.idx.nc(t,j)+1
0943
0944 mpc = mdi.flow(t,j,k).mpc;
0945 ion = find(mpc.branch(:, BR_STATUS));
0946 [Bdc, Bl, Psh, PLsh] = makeBdc(baseMVA, mpc.bus, mpc.branch(ion,:));
0947 mdi.flow(t,j,k).PLsh = PLsh;
0948 negCg = sparse(mpc.gen(:,GEN_BUS), (1:ng)', -1, ...
0949 mdi.idx.nb(t,j,k), ng);
0950 A = [Bdc negCg];
0951 b = -(mpc.bus(:,PD)+mpc.bus(:,GS))/baseMVA-Psh;
0952 vs = struct('name', {'Va', 'Pg'}, 'idx', {{t,j,k}, {t,j,k}});
0953 om.add_lin_constraint('Pmis', {t,j,k}, A, b, b, vs);
0954
0955 tmp = mpc.branch(ion,RATE_A)/baseMVA;
0956 iuncon = find(~tmp);
0957 tmp(iuncon) = Inf(size(iuncon));
0958 vs = struct('name', {'Va'}, 'idx', {{t,j,k}});
0959 om.add_lin_constraint('Pf', {t,j,k}, Bl, -tmp-PLsh, tmp-PLsh, vs);
0960 end
0961 end
0962 end
0963 else
0964 if verbose
0965 fprintf(' - Building load balance constraints.\n');
0966 end
0967
0968 for t = 1:nt
0969 for j = 1:mdi.idx.nj(t)
0970 for k = 1:mdi.idx.nc(t,j)+1
0971 mpc = mdi.flow(t,j,k).mpc;
0972 A = sparse(ones(1, ng));
0973 b = 1.0*sum(mpc.bus(:, PD)+mpc.bus(:,GS))/baseMVA;
0974 vs = struct('name', {'Pg'}, 'idx', {{t,j,k}});
0975 om.add_lin_constraint('Pmis', {t,j,k}, A, b, b, vs);
0976 end
0977 end
0978 end
0979 end
0980 if mdi.IncludeFixedReserves
0981 if verbose
0982 fprintf(' - Building fixed zonal reserve constraints.\n');
0983 end
0984 om.init_indexed_name('lin', 'Pg_plus_R', {nt, nj_max, nc_max+1});
0985 om.init_indexed_name('lin', 'Rreq', {nt, nj_max, nc_max+1});
0986 for t = 1:nt
0987 for j = 1:mdi.idx.nj(t)
0988 for k = 1:mdi.idx.nc(t,j)+1
0989
0990 mpc = mdi.flow(t,j,k).mpc;
0991 r = mdi.FixedReserves(t,j,k);
0992 ngr = length(r.igr);
0993 I = speye(ngr);
0994 Ar = sparse(1:ngr, r.igr, 1, ngr, ng);
0995 if UC
0996 A = [Ar I ...
0997 sparse(1:ngr, r.igr, -mpc.gen(r.igr, PMAX) / baseMVA, ngr, ng)];
0998 u = zeros(ngr, 1);
0999 vs = struct('name', {'Pg', 'R', 'u'}, 'idx', {{t,j,k}, {t,j,k}, {t}});
1000 else
1001 A = [Ar I];
1002 u = mpc.gen(r.igr, PMAX) / baseMVA;
1003 vs = struct('name', {'Pg', 'R'}, 'idx', {{t,j,k}, {t,j,k}});
1004 end
1005 om.add_lin_constraint('Pg_plus_R', {t,j,k}, A, [], u, vs);
1006 A = r.zones(:, r.igr);
1007 l = r.req / mpc.baseMVA;
1008 vs = struct('name', {'R'}, 'idx', {{t,j,k}});
1009 om.add_lin_constraint('Rreq', {t,j,k}, A, l, [], vs);
1010 end
1011 end
1012 end
1013 end
1014
1015
1016
1017 if verbose && ~isempty(mdi.Storage.UnitIdx)
1018 fprintf(' - Splitting storage injections into charge/discharge.\n');
1019 end
1020 om.init_indexed_name('lin', 'Ps', {nt, nj_max, nc_max+1});
1021 for t = 1:nt
1022 for j = 1:mdi.idx.nj(t)
1023 for k = 1:mdi.idx.nc(t,j)+1
1024 A = [sparse((1:ns)', mdi.Storage.UnitIdx, -1, ns, ng) Ins Ins];
1025 b = zeros(ns, 1);
1026 vs = struct('name', {'Pg', 'Psc', 'Psd'}, 'idx', {{t,j,k}, {t,j,k}, {t,j,k}});
1027 om.add_lin_constraint('Ps', {t,j,k}, A, b, b, vs);
1028 end
1029 end
1030 end
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040 if verbose
1041 fprintf(' - Building CCV constraints for piecewise-linear costs.\n');
1042 end
1043 om.init_indexed_name('lin', 'ycon', {nt, nj_max, nc_max+1});
1044 for t = 1:nt
1045 for j = 1:mdi.idx.nj(t)
1046 for k = 1:mdi.idx.nc(t,j)+1
1047 mpc = mdi.flow(t,j,k).mpc;
1048 [A, u] = makeAy(baseMVA, ng, mpc.gencost, 1, [], ng+1);
1049 vs = struct('name', {'Pg', 'y'}, 'idx', {{t,j,k}, {t,j,k}});
1050 om.add_lin_constraint('ycon', {t,j,k}, A, [], u, vs);
1051 end
1052 end
1053 end
1054
1055
1056
1057
1058
1059
1060
1061
1062 if verbose
1063 fprintf(' - Building contingency reserve constraints.\n');
1064 end
1065 om.init_indexed_name('lin', 'rampcont', {nt, nj_max, nc_max+1});
1066 for t =1:nt
1067 for j = 1:mdi.idx.nj(t)
1068 for k = 2:mdi.idx.nc(t,j)+1
1069 ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
1070 ngtmp = length(ii);
1071 A = sparse((1:ngtmp)', ii, 1, ngtmp, ng);
1072 u = mdi.flow(t,j,k).mpc.gen(ii,RAMP_10)/baseMVA;
1073 vs = struct('name', {'Pg', 'Pg'}, 'idx', {{t,j,1}, {t,j,k}});
1074 om.add_lin_constraint('rampcont', {t,j,k}, [-A A], -u, u, vs);
1075 end
1076 end
1077 end
1078
1079
1080
1081
1082
1083
1084
1085 om.init_indexed_name('lin', 'dPpRp', {nt, nj_max, nc_max+1});
1086 for t = 1:nt
1087 for j = 1:mdi.idx.nj(t);
1088 for k = 1:mdi.idx.nc(t,j)+1
1089 ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
1090 ngtmp = length(ii);
1091 A = sparse((1:ngtmp)', ii, 1, ngtmp, ng);
1092 l = zeros(ngtmp, 1);
1093 vs = struct('name', {'dPp', 'Rpp'}, 'idx', {{t,j,k}, {t}});
1094 om.add_lin_constraint('dPpRp', {t,j,k}, [-A A], l, [], vs);
1095 end
1096 end
1097 end
1098
1099
1100
1101
1102 om.init_indexed_name('lin', 'dPmRm', {nt, nj_max, nc_max+1});
1103 for t = 1:nt
1104 for j = 1:mdi.idx.nj(t);
1105 for k = 1:mdi.idx.nc(t,j)+1
1106 ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
1107 ngtmp = length(ii);
1108 A = sparse((1:ngtmp)', ii, 1, ngtmp, ng);
1109 l = zeros(ngtmp, 1);
1110 vs = struct('name', {'dPm', 'Rpm'}, 'idx', {{t,j,k}, {t}});
1111 om.add_lin_constraint('dPmRm', {t,j,k}, [-A A], l, [], vs);
1112 end
1113 end
1114 end
1115
1116
1117
1118 om.init_indexed_name('lin', 'dPdef', {nt, nj_max, nc_max+1});
1119 for t = 1:nt
1120 for j = 1:mdi.idx.nj(t)
1121 for k = 1:mdi.idx.nc(t,j)+1
1122 ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
1123 ngtmp = length(ii);
1124 A = sparse((1:ngtmp)', ii, 1, ngtmp, ng);
1125 b = zeros(ngtmp, 1);
1126 vs = struct('name', {'Pg', 'Pc', 'dPp', 'dPm'}, ...
1127 'idx', {{t,j,k}, {t}, {t,j,k}, {t,j,k}});
1128 om.add_lin_constraint('dPdef', {t,j,k}, [A -A -A A], b, b, vs);
1129 end
1130 end
1131 end
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145 if verbose
1146 fprintf(' - Building ramping transitions and reserve constraints.\n');
1147 end
1148 om.init_indexed_name('lin', 'Rrp', {nt, nj_max, nj_max});
1149
1150
1151 for t = 1:nt-1
1152 for j1 = 1:mdi.idx.nj(t)
1153 for j2 = 1:mdi.idx.nj(t+1)
1154 if mdi.tstep(t+1).TransMask(j2,j1)
1155 A = [Ing -Ing Ing];
1156 l = zeros(ng, 1);
1157 vs = struct('name', {'Pg', 'Pg', 'Rrp'}, ...
1158 'idx', {{t,j1,1}, {t+1,j2,1}, {t}});
1159 om.add_lin_constraint('Rrp', {t,j1,j2}, A, l, [], vs);
1160 end
1161 end
1162 end
1163 end
1164
1165
1166
1167
1168
1169 if ~mdi.OpenEnded
1170
1171 for j1 = 1:mdi.idx.nj(nt)
1172 A = [Ing Ing];
1173 l = mdi.TerminalPg/baseMVA;
1174 vs = struct('name', {'Pg', 'Rrp'}, ...
1175 'idx', {{nt,j1,1}, {nt}});
1176 om.add_lin_constraint('Rrp', {nt,j1,1}, A, l, [], vs);
1177 end
1178 end
1179
1180
1181
1182 om.init_indexed_name('lin', 'Rrm', {nt, nj_max, nj_max});
1183
1184
1185 for t = 1:nt-1
1186 for j1 = 1:mdi.idx.nj(t)
1187 for j2 = 1:mdi.idx.nj(t+1)
1188 if mdi.tstep(t+1).TransMask(j2,j1)
1189 A = [-Ing Ing Ing];
1190 l = zeros(ng, 1);
1191 vs = struct('name', {'Pg', 'Pg', 'Rrm'}, ...
1192 'idx', {{t,j1,1}, {t+1,j2,1}, {t}});
1193 om.add_lin_constraint('Rrm', {t,j1,j2}, A, l, [], vs);
1194 end
1195 end
1196 end
1197 end
1198
1199
1200
1201
1202
1203 if ~mdi.OpenEnded
1204
1205 for j1 = 1:mdi.idx.nj(nt)
1206 A = [-Ing Ing];
1207 l = -mdi.TerminalPg/baseMVA;
1208 vs = struct('name', {'Pg', 'Rrm'}, ...
1209 'idx', {{nt,j1,1}, {nt}});
1210 om.add_lin_constraint('Rrm', {nt,j1,1}, A, l, [], vs);
1211 end
1212 end
1213
1214
1215
1216 if ns
1217 if verbose
1218 fprintf(' - Building storage constraints.\n');
1219 end
1220
1221
1222 om.init_indexed_name('lin', 'Sm', {nt, nj_max});
1223 if mdi.Storage.ForceCyclicStorage
1224
1225 for j = 1:mdi.idx.nj(1)
1226 A = [ diagBeta2EtaIn1 diagBeta2overEtaOut1 Ins -spdiags(beta1(:,1), 0, ns, ns)];
1227 u = zeros(ns, 1);
1228 vs = struct('name', {'Psc', 'Psd', 'Sm', 'S0'}, 'idx', {{1,j,1}, {1,j,1}, {1}, {}});
1229 om.add_lin_constraint('Sm', {1,j}, A, [], u, vs);
1230 end
1231 else
1232
1233 for j = 1:mdi.idx.nj(1)
1234 A = [ diagBeta2EtaIn1 diagBeta2overEtaOut1 Ins ];
1235 u = beta1(:,1) .* ( rho(:,t).*mdi.Storage.InitialStorageLowerBound + ...
1236 (1-rho(:,t)).*mdi.Storage.InitialStorage ) / baseMVA;
1237 vs = struct('name', {'Psc', 'Psd', 'Sm'}, 'idx', {{1,j,1}, {1,j,1}, {1}});
1238 om.add_lin_constraint('Sm', {1,j}, A, [], u, vs);
1239 end
1240 end
1241
1242
1243
1244 for t = 2:nt
1245 for j = 1:mdi.idx.nj(t)
1246 Mj = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :) + ...
1247 mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1248 Lij = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1249 diag1minusRhoBeta1 = spdiags((1-rho(:,t)) .* beta1(:,t), 0, ns, ns);
1250 A = sparse([1:ns,1:ns,1:ns,1:ns]', ...
1251 [vv.i1.Psc(t,j,1):vv.iN.Psc(t,j,1), vv.i1.Psd(t,j,1):vv.iN.Psd(t,j,1), vv.i1.Sm(t-1):vv.iN.Sm(t-1), vv.i1.Sm(t):vv.iN.Sm(t)]', ...
1252 [beta2EtaIn(:,t); beta2overEtaOut(:,t); -beta1(:,t).*rho(:,t); ones(ns,1)], ...
1253 ns, nvars) ...
1254 - diag1minusRhoBeta1 * Mj;
1255 if mdi.Storage.ForceCyclicStorage
1256 As0 = sparse(ns, nvars);
1257 As0(:, vv.i1.S0:vv.iN.S0) = -diag1minusRhoBeta1 * Lij;
1258 A = A + As0;
1259 u = zeros(ns, 1);
1260 else
1261 u = full(diag1minusRhoBeta1 * Lij * mdi.Storage.InitialStorage/baseMVA);
1262 end
1263 om.add_lin_constraint('Sm', {t,j}, A, [], u);
1264 end
1265 end
1266
1267 om.init_indexed_name('lin', 'Sp', {nt, nj_max});
1268 if mdi.Storage.ForceCyclicStorage
1269
1270 for j = 1:mdi.idx.nj(1)
1271 A = [ -diagBeta2EtaIn1 -diagBeta2overEtaOut1 -Ins spdiags(beta1(:,1), 0, ns, ns) ];
1272 u = zeros(ns, 1);
1273 vs = struct('name', {'Psc', 'Psd', 'Sp', 'S0'}, 'idx', {{1,j,1}, {1,j,1}, {1}, {}});
1274 om.add_lin_constraint('Sp', {1,j}, A, [], u, vs);
1275 end
1276 else
1277
1278 for j = 1:mdi.idx.nj(1)
1279 A = [ -diagBeta2EtaIn1 -diagBeta2overEtaOut1 -Ins ];
1280 u = -beta1(:,1) .* ( rho(:,t).*mdi.Storage.InitialStorageUpperBound + ...
1281 (1-rho(:,t)).*mdi.Storage.InitialStorage ) / baseMVA;
1282 vs = struct('name', {'Psc', 'Psd', 'Sp'}, 'idx', {{1,j,1}, {1,j,1}, {1}});
1283 om.add_lin_constraint('Sp', {1,j}, A, [], u, vs);
1284 end
1285 end
1286
1287
1288
1289 for t = 2:nt
1290 for j = 1:mdi.idx.nj(t)
1291 Mj = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :) + ...
1292 mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1293 Lij = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1294 diag1minusRhoBeta1 = spdiags((1-rho(:,t)) .* beta1(:,t), 0, ns, ns);
1295 A = sparse([1:ns,1:ns,1:ns,1:ns]', ...
1296 [vv.i1.Psc(t,j,1):vv.iN.Psc(t,j,1), vv.i1.Psd(t,j,1):vv.iN.Psd(t,j,1), vv.i1.Sp(t-1):vv.iN.Sp(t-1), vv.i1.Sp(t):vv.iN.Sp(t)]', ...
1297 [-beta2EtaIn(:,t); -beta2overEtaOut(:,t); beta1(:,t).*rho(:,t); -ones(ns,1)], ...
1298 ns, nvars) ...
1299 + diag1minusRhoBeta1 * Mj;
1300 if mdi.Storage.ForceCyclicStorage
1301 As0 = sparse(ns, nvars);
1302 As0(:, vv.i1.S0:vv.iN.S0) = diag1minusRhoBeta1 * Lij;
1303 A = A + As0;
1304 u = zeros(ns, 1);
1305 else
1306 u = full(-diag1minusRhoBeta1 * Lij * mdi.Storage.InitialStorage/baseMVA);
1307 end
1308 om.add_lin_constraint('Sp', {t,j}, A, [], u);
1309 end
1310 end
1311
1312
1313
1314 om.init_indexed_name('lin', 'contSm', {nt, nj_max, nc_max+1});
1315 for j = 1:mdi.idx.nj(1)
1316 for k = 2:mdi.idx.nc(1,j)+1
1317 if mdi.Storage.ForceCyclicStorage
1318
1319 A = [ diagBeta4EtaIn1 diagBeta4overEtaOut1 diagBeta3EtaIn1 diagBeta3overEtaOut1 -spdiags(beta5(:,1), 0, ns, ns) ];
1320 u = -MinStorageLevel(:,1)/baseMVA;
1321 vs = struct('name', {'Psc', 'Psd', 'Psc', 'Psd', 'S0'}, ...
1322 'idx', {{1,j,1}, {1,j,1}, {1,j,k}, {1,j,k}, {}});
1323 else
1324
1325 A = [ diagBeta4EtaIn1 diagBeta4overEtaOut1 diagBeta3EtaIn1 diagBeta3overEtaOut1 ];
1326 u = (beta5(:,1).*mdi.Storage.InitialStorageLowerBound - MinStorageLevel(:,1))/baseMVA;
1327 vs = struct('name', {'Psc', 'Psd', 'Psc', 'Psd'}, ...
1328 'idx', {{1,j,1}, {1,j,1}, {1,j,k}, {1,j,k}});
1329 end
1330 om.add_lin_constraint('contSm', {1,j,k}, A, [], u, vs);
1331 end
1332 end
1333
1334
1335
1336 for t = 2:nt
1337 for j = 1:mdi.idx.nj(t)
1338 for k = 2:mdi.idx.nc(t,j)+1
1339 Mj = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :) + ...
1340 mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1341 Lij = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1342 diag1minusRhoBeta5 = spdiags((1-rho(:,t)) .* beta5(:,t), 0, ns, ns);
1343 A = sparse([1:ns,1:ns,1:ns,1:ns,1:ns]', ...
1344 [vv.i1.Psc(t,j,1):vv.iN.Psc(t,j,1), vv.i1.Psd(t,j,1):vv.iN.Psd(t,j,1), vv.i1.Psc(t,j,k):vv.iN.Psc(t,j,k), vv.i1.Psd(t,j,k):vv.iN.Psd(t,j,k), vv.i1.Sm(t-1):vv.iN.Sm(t-1)]', ...
1345 [beta4EtaIn(:,t); beta4overEtaOut(:,t); beta3EtaIn(:,t); beta3overEtaOut(:,t); -beta5(:,t).*rho(:,t)], ...
1346 ns, nvars) ...
1347 - diag1minusRhoBeta5 * Mj;
1348 u = -MinStorageLevel(:,t)/baseMVA;
1349 if mdi.Storage.ForceCyclicStorage
1350 As0 = sparse(ns, nvars);
1351 As0(:, vv.i1.S0:vv.iN.S0) = -diag1minusRhoBeta5 * Lij;
1352 A = A + As0;
1353 else
1354 u = u + diag1minusRhoBeta5 * Lij * mdi.Storage.InitialStorageLowerBound/baseMVA;
1355 end
1356 om.add_lin_constraint('contSm', {t,j,k}, A, [], u);
1357 end
1358 end
1359 end
1360
1361
1362 om.init_indexed_name('lin', 'contSp', {nt, nj_max, nc_max+1});
1363 for j = 1:mdi.idx.nj(1)
1364 for k = 2:mdi.idx.nc(1,j)+1
1365 if mdi.Storage.ForceCyclicStorage
1366
1367 A = [ -diagBeta4EtaIn1 -diagBeta4overEtaOut1 -diagBeta3EtaIn1 -diagBeta3overEtaOut1 spdiags(beta5(:,1), 0, ns, ns)];
1368 u = MaxStorageLevel(:,1)/baseMVA;
1369 vs = struct('name', {'Psc', 'Psd', 'Psc', 'Psd', 'S0'}, ...
1370 'idx', {{1,j,1}, {1,j,1}, {1,j,k}, {1,j,k}, {}});
1371 else
1372
1373 A = [ -diagBeta4EtaIn1 -diagBeta4overEtaOut1 -diagBeta3EtaIn1 -diagBeta3overEtaOut1 ];
1374 u = (MaxStorageLevel(:,1) - beta5(:,1).*mdi.Storage.InitialStorageUpperBound)/baseMVA;
1375 vs = struct('name', {'Psc', 'Psd', 'Psc', 'Psd'}, ...
1376 'idx', {{1,j,1}, {1,j,1}, {1,j,k}, {1,j,k}});
1377 end
1378 om.add_lin_constraint('contSp', {1,j,k}, A, [], u, vs);
1379 end
1380 end
1381
1382
1383
1384 for t = 2:nt
1385 for j = 1:mdi.idx.nj(t)
1386 for k = 2:mdi.idx.nc(t,j)+1
1387 Mj = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :) + ...
1388 mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1389 Lij = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1390 diag1minusRhoBeta5 = spdiags((1-rho(:,t)) .* beta5(:,t), 0, ns, ns);
1391 A = sparse([1:ns,1:ns,1:ns,1:ns,1:ns]', ...
1392 [vv.i1.Psc(t,j,1):vv.iN.Psc(t,j,1), vv.i1.Psd(t,j,1):vv.iN.Psd(t,j,1), vv.i1.Psc(t,j,k):vv.iN.Psc(t,j,k), vv.i1.Psd(t,j,k):vv.iN.Psd(t,j,k), vv.i1.Sp(t-1):vv.iN.Sp(t-1)]', ...
1393 [-beta4EtaIn(:,t); -beta4overEtaOut(:,t); -beta3EtaIn(:,t); -beta3overEtaOut(:,t); beta5(:,t).*rho(:,t)], ...
1394 ns, nvars) ...
1395 + diag1minusRhoBeta5 * Mj;
1396 u = MaxStorageLevel(:,t)/baseMVA;
1397 if mdi.Storage.ForceCyclicStorage
1398 As0 = sparse(ns, nvars);
1399 As0(:, vv.i1.S0:vv.iN.S0) = diag1minusRhoBeta5 * Lij;
1400 A = A + As0;
1401 else
1402 u = u - diag1minusRhoBeta5 * Lij * mdi.Storage.InitialStorageUpperBound/baseMVA;
1403 end
1404 om.add_lin_constraint('contSp', {t,j,k}, A, [], u);
1405 end
1406 end
1407 end
1408 end
1409
1410
1411
1412 if mdi.Storage.ForceExpectedTerminalStorage && mdi.Storage.ForceCyclicStorage
1413 error('most: ForceExpectedTerminalStorage and ForceCyclicStorage cannot be simultaneously true.');
1414 end
1415 if ns
1416
1417 if mdi.Storage.ForceExpectedTerminalStorage
1418
1419 A = sparse(ns, nvars);
1420 b = zeros(ns, 1);
1421 for j = 1:mdi.idx.nj(nt)
1422 Ngj = mdi.tstep(nt).Ng( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1423 Nhj = mdi.tstep(nt).Nh( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1424 Lfj = mdi.tstep(nt).Lf( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1425 A = A + mdi.CostWeights(1,j,nt) * (Ngj + Nhj);
1426 b = b + mdi.CostWeights(1,j,nt) * (Lfj * mdi.Storage.InitialStorage) / baseMVA;
1427 end
1428 endprob = sum(mdi.CostWeights(1,1:mdi.idx.nj(nt),nt)');
1429 A = (1/endprob) * A;
1430 b = (1/endprob) * b;
1431 l = mdi.Storage.ExpectedTerminalStorageMin / baseMVA - b;
1432 u = mdi.Storage.ExpectedTerminalStorageMax / baseMVA - b;
1433 om.add_lin_constraint('ESnt', A, l, u);
1434 elseif mdi.Storage.ForceCyclicStorage
1435
1436 A = sparse(ns, nvars);
1437 for j = 1:mdi.idx.nj(nt)
1438 Ngj = mdi.tstep(nt).Ng( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1439 Nhj = mdi.tstep(nt).Nh( j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1440 A = A + mdi.CostWeights(1,j,nt) * (Ngj + Nhj);
1441 end
1442 endprob = sum(mdi.CostWeights(1,1:mdi.idx.nj(nt),nt)');
1443 A = (1/endprob) * A;
1444 for j = 1:mdi.idx.nj(nt)
1445 Lfj = mdi.tstep(nt).Lf(j:mdi.idx.nj(nt):(ns-1)*mdi.idx.nj(nt)+j, :);
1446 A(:, vv.i1.S0:vv.iN.S0) = A(:, vv.i1.S0:vv.iN.S0) ...
1447 + mdi.CostWeights(1,j,nt) * Lfj;
1448 end
1449 A(:, vv.i1.S0:vv.iN.S0) = (1/endprob) * A(:, vv.i1.S0:vv.iN.S0) - speye(ns);
1450 b = zeros(ns, 1);
1451 om.add_lin_constraint('ESnt', A, b, b);
1452 end
1453 end
1454
1455
1456 if nzds || nyds
1457 if verbose
1458 fprintf(' - Building dynamical system constraints.\n');
1459 end
1460
1461
1462
1463 for t = 1:nt
1464 E = sparse(ng, nvars);
1465 for j = 1:mdi.idx.nj(t)
1466 for k = 1:mdi.idx.nc(t,j)+1;
1467 E = E + (mdi.CostWeightsAdj(k,j,t)/mdi.StepProb(t)) * sparse((1:ng)', (vv.i1.Pg(t,j,k):vv.iN.Pg(t,j,k))', 1, ng, nvars);
1468 end
1469 end
1470 mdi.tstep(t).E = E;
1471 end
1472 end
1473
1474
1475
1476 if nzds
1477 om.init_indexed_name('lin', 'DSz', {ntds-1});
1478 b = zeros(nzds, 1);
1479 for t = 1:ntds-1
1480 if t <= nt
1481
1482 A = mdi.dstep(t).B * mdi.tstep(t).E;
1483 else
1484
1485
1486
1487
1488 A = sparse(nzds, nvars);
1489 end
1490 A(:, vv.i1.Z(t):vv.iN.Z(t)) = mdi.dstep(t).A;
1491 A(:, vv.i1.Z(t+1):vv.iN.Z(t+1)) = -speye(nzds);
1492 om.add_lin_constraint('DSz', {t}, A, b, b);
1493 end
1494 end
1495
1496
1497 if nyds
1498 om.init_indexed_name('lin', 'DSy', {ntds});
1499 for t = 1:ntds
1500 if t <= nt
1501 A = mdi.dstep(t).D * mdi.tstep(t).E;
1502 else
1503 A = sparse(nyds, nvars);
1504 end
1505 if nzds
1506 A(:, vv.i1.Z(t):vv.iN.Z(t)) = mdi.dstep(t).C;
1507 end
1508 l = mdi.dstep(t).ymin;
1509 u = mdi.dstep(t).ymax;
1510 om.add_lin_constraint('DSy', {t}, A, l, u);
1511 end
1512 end
1513
1514
1515 if UC
1516 if verbose
1517 fprintf(' - Building unit commitment constraints.\n');
1518 end
1519
1520 om.init_indexed_name('lin', 'uvw', {nt});
1521 for t = 1:nt
1522 if t == 1
1523
1524 if mdi.UC.CyclicCommitment
1525 vs = struct('name', {'u', 'u', 'v', 'w'}, 'idx', {{1}, {nt}, {1}, {1}});
1526 A = [Ing -Ing -Ing Ing];
1527 b = zeros(ng, 1);
1528 else
1529 vs = struct('name', {'u', 'v', 'w'}, 'idx', {{1}, {1}, {1}});
1530 A = [Ing -Ing Ing];
1531 b = (mdi.UC.InitialState > 0);
1532 end
1533 else
1534
1535 vs = struct('name', {'u', 'u', 'v', 'w'}, 'idx', {{t}, {t-1}, {t}, {t}});
1536 A = [Ing -Ing -Ing Ing];
1537 b = zeros(ng, 1);
1538 end
1539 om.add_lin_constraint('uvw', {t}, A, b, b, vs);
1540 end
1541
1542
1543 om.init_indexed_name('lin', 'minup', {nt, ng});
1544 for t = 1:nt
1545 for i = 1:ng
1546 ti = t-mdi.UC.MinUp(i)+1:t;
1547 if mdi.UC.CyclicCommitment
1548 for tt = 1:length(ti)
1549 if ti(tt) < 1
1550 ti(tt) = nt + ti(tt);
1551 end
1552 end
1553 end
1554
1555
1556
1557 ti = ti(ti>0);
1558 vs = struct('name', {'u'}, 'idx', {{t}});
1559 A = sparse(1, i, -1, 1, ng);
1560 for tt = 1:length(ti)
1561 vs(end+1).name = 'v';
1562 vs(end).idx = {ti(tt)};
1563 A = [A sparse(1, i, 1, 1, ng)];
1564 end
1565 om.add_lin_constraint('minup', {t, i}, A, [], 0, vs);
1566 end
1567 end
1568
1569
1570 om.init_indexed_name('lin', 'mindown', {nt, ng});
1571 for t = 1:nt
1572 for i = 1:ng
1573 ti = t-mdi.UC.MinDown(i)+1:t;
1574 if mdi.UC.CyclicCommitment
1575 for tt = 1:length(ti)
1576 if ti(tt) < 1
1577 ti(tt) = nt + ti(tt);
1578 end
1579 end
1580 end
1581
1582
1583
1584 ti = ti(ti>0);
1585 vs = struct('name', {'u'}, 'idx', {{t}});
1586 A = sparse(1, i, 1, 1, ng);
1587 for tt = 1:length(ti)
1588 vs(end+1).name = 'w';
1589 vs(end).idx = {ti(tt)};
1590 A = [A sparse(1, i, 1, 1, ng)];
1591 end
1592 om.add_lin_constraint('mindown', {t, i}, A, [], 1, vs);
1593 end
1594 end
1595
1596
1597
1598
1599 om.init_indexed_name('lin', 'uPmax', {nt, nj_max, nc_max+1});
1600 for t = 1:nt
1601 for j = 1:mdi.idx.nj(t)
1602 for k = 1:mdi.idx.nc(t,j)+1
1603 mpc = mdi.flow(t,j,k).mpc;
1604 ii = find(mpc.gen(:, GEN_STATUS));
1605 nii = length(ii);
1606 vs = struct('name', {'Pg', 'u'}, 'idx', {{t,j,k}, {t}});
1607 A = [ sparse(1:nii, ii, 1, nii, ng) ...
1608 sparse(1:nii, ii, -mpc.gen(ii, PMAX)/baseMVA, nii, ng) ];
1609 u = zeros(nii, 1);
1610 om.add_lin_constraint('uPmax', {t,j,k}, A, [], u, vs);
1611 end
1612 end
1613 end
1614
1615 om.init_indexed_name('lin', 'uPmin', {nt, nj_max, nc_max+1});
1616 for t = 1:nt
1617 for j = 1:mdi.idx.nj(t)
1618 for k = 1:mdi.idx.nc(t,j)+1
1619 mpc = mdi.flow(t,j,k).mpc;
1620 ii = find(mpc.gen(:, GEN_STATUS));
1621 nii = length(ii);
1622 vs = struct('name', {'Pg', 'u'}, 'idx', {{t,j,k}, {t}});
1623 A = [ sparse(1:nii, ii, -1, nii, ng) ...
1624 sparse(1:nii, ii, mpc.gen(ii, PMIN)/baseMVA, nii, ng) ];
1625 u = zeros(nii, 1);
1626 om.add_lin_constraint('uPmin', {t,j,k}, A, [], u, vs);
1627 end
1628 end
1629 end
1630
1631
1632
1633
1634 if mdi.QCoordination
1635 om.init_indexed_name('lin', 'uQmax', {nt, nj_max, nc_max+1});
1636 for t = 1:nt
1637 for j = 1:mdi.idx.nj(t)
1638 for k = 1:mdi.idx.nc(t,j)+1
1639 mpc = mdi.flow(t,j,k).mpc;
1640 ii = find(mpc.gen(:, GEN_STATUS));
1641 nii = length(ii);
1642 vs = struct('name', {'Qg', 'u'}, 'idx', {{t,j,k}, {t}});
1643 A = [ sparse(1:nii, ii, 1, nii, ng) ...
1644 sparse(1:nii, ii, -mpc.gen(ii, QMAX)/baseMVA, nii, ng) ];
1645 u = zeros(nii, 1);
1646 om.add_lin_constraint('uQmax', {t,j,k}, A, [], u, vs);
1647 end
1648 end
1649 end
1650
1651 om.init_indexed_name('lin', 'uQmin', {nt, nj_max, nc_max+1});
1652 for t = 1:nt
1653 for j = 1:mdi.idx.nj(t)
1654 for k = 1:mdi.idx.nc(t,j)+1
1655 mpc = mdi.flow(t,j,k).mpc;
1656 ii = find(mpc.gen(:, GEN_STATUS));
1657 nii = length(ii);
1658 vs = struct('name', {'Qg', 'u'}, 'idx', {{t,j,k}, {t}});
1659 A = [ sparse(1:nii, ii, -1, nii, ng) ...
1660 sparse(1:nii, ii, mpc.gen(ii, QMIN)/baseMVA, nii, ng) ];
1661 u = zeros(nii, 1);
1662 om.add_lin_constraint('uQmin', {t,j,k}, A, [], u, vs);
1663 end
1664 end
1665 end
1666 end
1667 end
1668
1669 if verbose
1670 fprintf('- Building cost structures.\n');
1671 end
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685 if mdi.OpenEnded
1686 om.init_indexed_name('qdc', 'RampWear', {nt, nj_max, nj_max});
1687 else
1688 om.init_indexed_name('qdc', 'RampWear', {nt+1, nj_max, nj_max});
1689 end
1690 for j = 1:mdi.idx.nj(1)
1691 w = mdi.tstep(1).TransMat(j,1);
1692 Q = spdiags(w * baseMVA^2 * mdi.RampWearCostCoeff(:,1), 0, ng, ng);
1693 c = -w * baseMVA * mdi.RampWearCostCoeff(:,1) .* mdi.InitialPg;
1694 vs = struct('name', {'Pg'}, 'idx', {{1,j,1}});
1695 k0 = w * 0.5 * mdi.RampWearCostCoeff(:,1)' * mdi.InitialPg.^2;
1696 om.add_quad_cost('RampWear', {1,j,1}, Q, c, k0, vs);
1697 end
1698
1699 for t = 2:nt
1700 for j2 = 1:mdi.idx.nj(t)
1701 for j1 = 1:mdi.idx.nj(t-1)
1702 w = mdi.tstep(t).TransMat(j2,j1) * mdi.CostWeights(1, j1, t-1);
1703 h = w * baseMVA^2 * mdi.RampWearCostCoeff(:,t);
1704 i = (1:ng)';
1705 j = ng+(1:ng)';
1706 Q = sparse([i;j;i;j], [i;i;j;j], [h;-h;-h;h], 2*ng, 2*ng);
1707 vs = struct('name', {'Pg', 'Pg'}, 'idx', {{t-1,j1,1}, {t,j2,1}});
1708 om.add_quad_cost('RampWear', {t,j1,j2}, Q, zeros(2*ng,1), 0, vs);
1709 end
1710 end
1711 end
1712
1713
1714
1715
1716 if ~mdi.OpenEnded
1717 for j = 1:mdi.idx.nj(nt)
1718 w = mdi.tstep(nt+1).TransMat(1, j) * mdi.CostWeights(1, j, nt);
1719 Q = spdiags(w * baseMVA^2 * mdi.RampWearCostCoeff(:,nt+1), 0, ng, ng);
1720 c = -w * baseMVA * mdi.RampWearCostCoeff(:,nt+1) .* mdi.TerminalPg;
1721 vs = struct('name', {'Pg'}, 'idx', {{nt,j,1}});
1722 k0 = w * 0.5 * mdi.RampWearCostCoeff(:,nt+1)' * mdi.TerminalPg.^2;
1723 om.add_quad_cost('RampWear', {nt+1,j,1}, Q, c, k0, vs);
1724 end
1725 end
1726
1727
1728
1729 om.init_indexed_name('qdc', 'Cp', {nt, nj_max, nc_max+1});
1730 om.init_indexed_name('qdc', 'Cy', {nt, nj_max, nc_max+1});
1731 om.init_indexed_name('qdc', 'Cpp', {nt, nj_max, nc_max+1});
1732 om.init_indexed_name('qdc', 'Cpm', {nt, nj_max, nc_max+1});
1733 if mdi.IncludeFixedReserves
1734 om.init_indexed_name('qdc', 'Rcost', {nt, nj_max, nc_max+1});
1735 end
1736 om.init_indexed_name('qdc', 'Crpp', {nt});
1737 om.init_indexed_name('qdc', 'Crpm', {nt});
1738 for t = 1:nt
1739 for j = 1:mdi.idx.nj(t)
1740 for k = 1:mdi.idx.nc(t,j)+1
1741 w = mdi.CostWeightsAdj(k,j,t);
1742
1743
1744 gc = mdi.flow(t,j,k).mpc.gencost;
1745 ipol = find(gc(:, MODEL) == POLYNOMIAL);
1746 if ~isempty(ipol)
1747 ncost = gc(ipol(1), NCOST);
1748 if all(gc(ipol, NCOST) == ncost)
1749
1750 if ncost > 3
1751 error('most: polynomial generator costs of order higher than quadratic not supported');
1752 elseif ncost == 3
1753 Q = sparse(ipol, ipol, 2 * w * baseMVA^2*gc(ipol, COST), ng, ng);
1754 else
1755 Q = sparse(ng,ng);
1756 end
1757 c = zeros(ng, 1);
1758 if ncost >= 2
1759 c(ipol) = w * baseMVA*gc(ipol, COST+ncost-2);
1760 end
1761 k0 = w * sum(gc(ipol, COST+ncost-1));
1762 else
1763
1764 Q = sparse(ng,ng);
1765 c = zeros(ng, 1);
1766 for i = ipol'
1767 ncost = gc(i, NCOST);
1768 if ncost > 3
1769 error('most: polynomial generator costs of order higher than quadratic not supported');
1770 elseif ncost == 3
1771 Q(i,i) = 2 * w * baseMVA^2*gc(i, COST);
1772 end
1773 if ncost >= 2
1774 c(i) = w * baseMVA*gc(i, COST+ncost-2);
1775 end
1776 k0 = w * gc(i, COST+ncost-1);
1777 end
1778 end
1779 vs = struct('name', {'Pg'}, 'idx', {{t,j,k}});
1780 om.add_quad_cost('Cp', {t,j,k}, Q, c, k0, vs);
1781 end
1782
1783
1784
1785 if mdi.idx.ny(t,j,k)
1786 c = w * ones(mdi.idx.ny(t,j,k),1);
1787 vs = struct('name', {'y'}, 'idx', {{t,j,k}});
1788 om.add_quad_cost('Cy', {t,j,k}, [], c, 0, vs);
1789 end
1790
1791
1792 c = w * baseMVA * mdi.offer(t).PositiveActiveDeltaPrice(:);
1793 vs = struct('name', {'dPp'}, 'idx', {{t,j,k}});
1794 om.add_quad_cost('Cpp', {t,j,k}, [], c, 0, vs);
1795 c = w * baseMVA * mdi.offer(t).NegativeActiveDeltaPrice(:);
1796 vs = struct('name', {'dPm'}, 'idx', {{t,j,k}});
1797 om.add_quad_cost('Cpm', {t,j,k}, [], c, 0, vs);
1798
1799
1800 if mdi.IncludeFixedReserves
1801 c = w * mdi.FixedReserves(t,j,k).cost(r.igr) * baseMVA;
1802 vs = struct('name', {'R'}, 'idx', {{t,j,k}});
1803 om.add_quad_cost('Rcost', {t,j,k}, [], c, 0, vs);
1804 end
1805 end
1806 end
1807
1808
1809 c = baseMVA * mdi.StepProb(t) * mdi.offer(t).PositiveActiveReservePrice(:);
1810 vs = struct('name', {'Rpp'}, 'idx', {{t}});
1811 om.add_quad_cost('Crpp', {t}, [], c, 0, vs);
1812 c = baseMVA * mdi.StepProb(t) * mdi.offer(t).NegativeActiveReservePrice(:);
1813 vs = struct('name', {'Rpm'}, 'idx', {{t}});
1814 om.add_quad_cost('Crpm', {t}, [], c, 0, vs);
1815 end
1816
1817 om.init_indexed_name('qdc', 'Crrp', {mdi.idx.ntramp});
1818 om.init_indexed_name('qdc', 'Crrm', {mdi.idx.ntramp});
1819 for t = 1:nt-1,
1820 c = baseMVA * mdi.StepProb(t+1) * mdi.offer(t).PositiveLoadFollowReservePrice(:);
1821 vs = struct('name', {'Rrp'}, 'idx', {{t}});
1822 om.add_quad_cost('Crrp', {t}, [], c, 0, vs);
1823 c = baseMVA * mdi.StepProb(t+1) * mdi.offer(t).NegativeLoadFollowReservePrice(:);
1824 vs = struct('name', {'Rrm'}, 'idx', {{t}});
1825 om.add_quad_cost('Crrm', {t}, [], c, 0, vs);
1826 end
1827
1828 if ~mdi.OpenEnded
1829
1830 c = baseMVA * mdi.offer(nt).PositiveLoadFollowReservePrice(:);
1831 vs = struct('name', {'Rrp'}, 'idx', {{nt}});
1832 om.add_quad_cost('Crrp', {nt}, [], c, 0, vs);
1833 c = baseMVA * mdi.offer(nt).NegativeLoadFollowReservePrice(:);
1834 vs = struct('name', {'Rrm'}, 'idx', {{nt}});
1835 om.add_quad_cost('Crrm', {nt}, [], c, 0, vs);
1836 end
1837
1838 if UC
1839 om.init_indexed_name('qdc', 'c00', {nt});
1840 om.init_indexed_name('qdc', 'startup', {nt});
1841 om.init_indexed_name('qdc', 'shutdown', {nt});
1842 for t = 1:nt
1843 ww = zeros(ng, 1);
1844 for j = 1:mdi.idx.nj(t)
1845 for k = 1:mdi.idx.nc(t)+1
1846 ww = ww + mdi.CostWeightsAdj(k,j,t) * mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS);
1847 end
1848 end
1849 c = ww.*mdi.UC.c00(:,t);
1850 vs = struct('name', {'u'}, 'idx', {{t}});
1851 om.add_quad_cost('c00', {t}, [], c, 0, vs);
1852 c = mdi.StepProb(t)*mdi.flow(t,1,1).mpc.gencost(:, STARTUP);
1853 vs = struct('name', {'v'}, 'idx', {{t}});
1854 om.add_quad_cost('startup', {t}, [], c, 0, vs);
1855 c = mdi.StepProb(t)*mdi.flow(t,1,1).mpc.gencost(:, SHUTDOWN);
1856 vs = struct('name', {'w'}, 'idx', {{t}});
1857 om.add_quad_cost('shutdown', {t}, [], c, 0, vs);
1858 end
1859 end
1860
1861 if ns
1862 A1 = sparse(ns, ns);
1863 A2 = sparse(ns, nvars);
1864 A3 = sparse(ns, nvars);
1865 A4 = sparse(ns, nvars);
1866 A5 = sparse(ns, nvars);
1867 A6 = sparse(ns, nvars);
1868 A7 = sparse(ns, nvars);
1869
1870 vv = om.get_idx();
1871 for t = 1:nt
1872
1873
1874 for j = 1:mdi.idx.nj(t)
1875
1876 Gtj0 = mdi.tstep(t).G( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1877 Htj0 = mdi.tstep(t).H( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1878 Litj0 = mdi.tstep(t).Li( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1879 Mgtj0 = mdi.tstep(t).Mg( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1880 Mhtj0 = mdi.tstep(t).Mh( j:mdi.idx.nj(t):(ns-1)*mdi.idx.nj(t)+j, :);
1881 sum_psi_tjk = sum(mdi.CostWeights(2:mdi.idx.nc(t,j)+1,j,t));
1882 if t == nt
1883 A1 = A1 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta1(:,t), 0, ns, ns) * Litj0;
1884 A2 = A2 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta1(:,t), 0, ns, ns) * Mgtj0;
1885 A3 = A3 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta1(:,t), 0, ns, ns) * Mhtj0;
1886 A4 = A4 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta2(:,t), 0, ns, ns) * Gtj0;
1887 A5 = A5 + mdi.CostWeights(1,j,t) * spdiags(OutEff(:,t) .* beta2(:,t), 0, ns, ns) * Htj0;
1888 end
1889 A1 = A1 + sum_psi_tjk * spdiags(OutEff(:,t) .* beta5(:,t), 0, ns, ns) * Litj0;
1890 A2 = A2 + sum_psi_tjk * (spdiags(OutEff(:,t) .* beta5(:,t), 0, ns, ns) * Mgtj0 + spdiags(OutEff(:,t) .* beta4(:,t), 0, ns, ns) * Gtj0);
1891 A3 = A3 + sum_psi_tjk * (spdiags(OutEff(:,t) .* beta5(:,t), 0, ns, ns) * Mhtj0 + spdiags(OutEff(:,t) .* beta4(:,t), 0, ns, ns) * Htj0);
1892 for k = 2:mdi.idx.nc(t,j)+1
1893 ii = (1:ns)';
1894 jj1 = (vv.i1.Psc(t,j,k):vv.iN.Psc(t,j,k))';
1895 jj2 = (vv.i1.Psd(t,j,k):vv.iN.Psd(t,j,k))';
1896 Gtjk = sparse(ii, jj1, -mdi.Delta_T * InEff(:,t), ns, nvars);
1897 Htjk = sparse(ii, jj2, -mdi.Delta_T ./ OutEff(:,t), ns, nvars);
1898 A6 = A6 + mdi.CostWeights(k,j,t) * spdiags(OutEff(:,t) .* beta3(:,t), 0, ns, ns) * Gtjk;
1899 A7 = A7 + mdi.CostWeights(k,j,t) * spdiags(OutEff(:,t) .* beta3(:,t), 0, ns, ns) * Htjk;
1900 end
1901 end
1902 end
1903 Cfstor = -baseMVA * ...
1904 (mdi.Storage.TerminalStoragePrice' * (A2 + A3) + ...
1905 mdi.Storage.TerminalChargingPrice0' * A4 + ...
1906 mdi.Storage.TerminalDischargingPrice0' * A5 + ...
1907 mdi.Storage.TerminalChargingPriceK' * A6 + ...
1908 mdi.Storage.TerminalDischargingPriceK' * A7);
1909 if mdi.Storage.ForceCyclicStorage
1910
1911
1912
1913
1914 Cfstor(vv.i1.S0:vv.iN.S0) = ...
1915 Cfstor(vv.i1.S0:vv.iN.S0) + ...
1916 baseMVA * mdi.Storage.InitialStorageCost';
1917
1918
1919 Cfstor(vv.i1.S0:vv.iN.S0) = ...
1920 Cfstor(vv.i1.S0:vv.iN.S0) - ...
1921 baseMVA * mdi.Storage.TerminalStoragePrice' * A1;
1922 end
1923 om.add_quad_cost('fstor', [], Cfstor', 0);
1924
1925
1926
1927 om.init_indexed_name('qdc', 'SpSmFudge', {nt});
1928 c = 1e-2 * [-ones(ns,1); ones(ns,1)];
1929 for t = 1:nt
1930 vs = struct('name', {'Sm', 'Sp'}, 'idx', {{t}, {t}});
1931 om.add_quad_cost('SpSmFudge', {t}, [], c, 0, vs);
1932 end
1933 else
1934 Cfstor = sparse(1, nvars);
1935 end
1936
1937
1938
1939
1940 if verbose
1941 fprintf('- Assembling full set of constraints.\n');
1942 end
1943 [mdi.QP.A, mdi.QP.l, mdi.QP.u] = om.params_lin_constraint();
1944 if verbose
1945 fprintf('- Assembling full set of variable bounds.\n');
1946 end
1947 [mdi.QP.x0, mdi.QP.xmin, mdi.QP.xmax, mdi.QP.vtype] = om.params_var();
1948 if verbose
1949 fprintf('- Assembling full set of costs.\n');
1950 end
1951 [mdi.QP.H1, mdi.QP.C1, mdi.QP.c1] = om.params_quad_cost();
1952
1953
1954 mdi.QP.Cfstor = Cfstor;
1955 mdi.om = om;
1956 else
1957 om = mdi.om;
1958 end
1959
1960
1961
1962 if isfield(mdi, 'CoordCost') && ...
1963 (~isempty(mdi.CoordCost.Cuser) || ~isempty(mdi.CoordCost.Huser))
1964 if verbose
1965 fprintf('- Adding coordination cost to standard cost.\n');
1966 end
1967 nvuser = length(mdi.CoordCost.Cuser);
1968 nvdiff = mdi.idx.nvars - nvuser;
1969 om.add_quad_cost( 'CoordCost', ...
1970 [ mdi.CoordCost.Huser sparse(nvuser,nvdiff);
1971 sparse(nvdiff,nvuser) sparse(nvdiff,nvdiff) ], ...
1972 mdi.CoordCost.Cuser(:), mdi.CoordCost.cuser);
1973 end
1974 [mdi.QP.H, mdi.QP.C, mdi.QP.c] = om.params_quad_cost();
1975
1976 et_setup = toc(t0);
1977 t0 = tic;
1978
1979
1980 mdo = mdi;
1981 mdo.om = om.copy();
1982
1983 if mpopt.most.solve_model
1984
1985 if mdi.DCMODEL ~= mo.DCMODEL
1986 error('MDI.DCMODEL inconsistent with MPOPT.most.dc_model');
1987 end
1988 if mdi.IncludeFixedReserves ~= mo.IncludeFixedReserves
1989 error('MDI.IncludeFixedReserves inconsistent with MPOPT.most.fixed_res (and possible presence of MDI.FixedReserves(t,j,k))');
1990 end
1991 if mdi.SecurityConstrained ~= mo.SecurityConstrained
1992 error('MDI.SecurityConstrained inconsistent with MPOPT.most.security_constraints (and possible presence of MDI.cont(t,j).contab)');
1993 end
1994 if mdi.QCoordination ~= mo.QCoordination
1995 error('MDI.QCoordination inconsistent with MPOPT.most.q_coordination');
1996 end
1997 if mdi.Storage.ForceCyclicStorage ~= mo.ForceCyclicStorage
1998 error('MDI.Storage.ForceCyclicStorage inconsistent with MPOPT.most.storage.cyclic');
1999 end
2000 if mdi.Storage.ForceExpectedTerminalStorage ~= mo.ForceExpectedTerminalStorage
2001 error('MDI.Storage.ForceExpectedTerminalStorage inconsistent with MPOPT.most.storage.terminal_target (and possible presence of MDI.Storage.ExpectedTerminalStorageAim|Min|Max)');
2002 end
2003 if mdi.UC.run ~= UC
2004 error('MDI.UC.run inconsistent with MPOPT.most.uc.run (and possible presence of MDI.UC.CommitKey)');
2005 end
2006
2007 model = om.problem_type();
2008 mdo.QP.opt = mpopt2qpopt(mpopt, model, 'most');
2009 mdo.QP.opt.x0 = [];
2010 if verbose
2011 fprintf('- Calling %s solver.\n\n', model);
2012 fprintf('============================================================================\n\n');
2013 end
2014 [mdo.QP.x, mdo.QP.f, mdo.QP.exitflag, mdo.QP.output, mdo.QP.lambda ] = ...
2015 mdo.om.solve(mdo.QP.opt);
2016 if mdo.QP.exitflag > 0
2017 success = 1;
2018 if verbose
2019 fprintf('\n============================================================================\n');
2020 fprintf('- MOST: %s solved successfully.\n', model);
2021 end
2022 else
2023 success = 0;
2024 if verbose
2025 fprintf('\n============================================================================\n');
2026 fprintf('- MOST: %s solver ''%s'' failed with exit flag = %d\n', model, mdo.QP.opt.alg, mdo.QP.exitflag);
2027 end
2028
2029
2030
2031
2032
2033 end
2034
2035 if verbose
2036 fprintf('- Post-processing results.\n');
2037 end
2038 if success
2039 om = mdo.om;
2040 for t = 1:nt
2041 if UC
2042 mdo.UC.CommitSched(:, t) = om.get_soln('var', 'u', {t});
2043 end
2044 for j = 1:mdi.idx.nj(t)
2045 for k = 1:mdi.idx.nc(t,j)+1
2046 mpc = mdo.flow(t,j,k).mpc;
2047
2048 if mdo.DCMODEL
2049 mpc.bus(:, VM) = 1;
2050 end
2051
2052 mpc.gen(:, PG) = baseMVA * om.get_soln('var', 'Pg', {t,j,k});
2053
2054 Pmin = mpc.gen(:, PMIN);
2055 Qmin = mpc.gen(:, QMIN);
2056 Qmax = mpc.gen(:, QMAX);
2057 ivl = find( isload(mpc.gen) & (Qmin ~= 0 | Qmax ~= 0) );
2058 Qlim = (Qmin(ivl) == 0) .* Qmax(ivl) + (Qmax(ivl) == 0) .* Qmin(ivl);
2059 mpc.gen(ivl, QG) = mpc.gen(ivl, PG) .* Qlim ./ Pmin(ivl);
2060 if mdo.DCMODEL
2061
2062 Va = om.get_soln('var', 'Va', {t,j,k});
2063 mpc.bus(:, VA) = (180/pi) * Va;
2064
2065
2066 price = (om.get_soln('lin', 'mu_u', 'Pmis', {t,j,k})-om.get_soln('lin', 'mu_l', 'Pmis', {t,j,k})) / baseMVA;
2067 mpc.bus(:, LAM_P) = price;
2068
2069
2070 mpc.branch(:, PF) = 0;
2071 mpc.branch(:, QF) = 0;
2072 mpc.branch(:, PT) = 0;
2073 mpc.branch(:, QT) = 0;
2074 mpc.branch(:, MU_SF) = 0;
2075 mpc.branch(:, MU_ST) = 0;
2076 ion = find(mpc.branch(:, BR_STATUS));
2077 AA = om.params_lin_constraint('Pf', {t,j,k});
2078 lf = baseMVA * (AA * Va + mdo.flow(t,j,k).PLsh);
2079 mpc.branch(ion, PF) = lf;
2080 mpc.branch(ion, PT) = -lf;
2081 mpc.branch(ion, MU_SF) = om.get_soln('lin', 'mu_u', 'Pf', {t,j,k}) / baseMVA;
2082 mpc.branch(ion, MU_ST) = om.get_soln('lin', 'mu_l', 'Pf', {t,j,k}) / baseMVA;
2083 else
2084
2085 price = (om.get_soln('lin', 'mu_l', 'Pmis', {t,j,k})-om.get_soln('lin', 'mu_u', 'Pmis', {t,j,k})) / baseMVA;
2086 mpc.bus(:, LAM_P) = price;
2087 end
2088 if UC
2089
2090
2091 igenon = find(mpc.gen(:, GEN_STATUS));
2092 mpc.gen(igenon, GEN_STATUS) = mdo.UC.CommitSched(igenon, t);
2093 gs = mpc.gen(igenon, GEN_STATUS) > 0;
2094 mpc.gen(:, MU_PMAX) = 0;
2095 mpc.gen(:, MU_PMIN) = 0;
2096 mpc.gen(igenon, MU_PMAX) = gs .* ...
2097 om.get_soln('lin', 'mu_u', 'uPmax', {t,j,k}) / baseMVA;
2098 mpc.gen(igenon, MU_PMIN) = gs .* ...
2099 om.get_soln('lin', 'mu_u', 'uPmin', {t,j,k}) / baseMVA;
2100 if mdo.QCoordination
2101 mpc.gen(:, MU_QMAX) = 0;
2102 mpc.gen(:, MU_QMIN) = 0;
2103 mpc.gen(igenon, MU_QMAX) = gs .* ...
2104 om.get_soln('lin', 'mu_u', 'uQmax', {t,j,k}) / baseMVA;
2105 mpc.gen(igenon, MU_QMIN) = gs .* ...
2106 om.get_soln('lin', 'mu_u', 'uQmin', {t,j,k}) / baseMVA;
2107 end
2108 else
2109 gs = mpc.gen(:, GEN_STATUS) > 0;
2110 mpc.gen(:, MU_PMAX) = gs .* ...
2111 om.get_soln('var', 'mu_u', 'Pg', {t,j,k}) / baseMVA;
2112 mpc.gen(:, MU_PMIN) = gs .* ...
2113 om.get_soln('var', 'mu_l', 'Pg', {t,j,k}) / baseMVA;
2114 if mdo.QCoordination
2115 mpc.gen(:, MU_QMAX) = gs .* ...
2116 om.get_soln('var', 'mu_u', 'Qg', {t,j,k}) / baseMVA;
2117 mpc.gen(:, MU_QMIN) = gs .* ...
2118 om.get_soln('var', 'mu_l', 'Qg', {t,j,k}) / baseMVA;
2119 end
2120 end
2121 if mdi.IncludeFixedReserves
2122 z = zeros(ng, 1);
2123 r = mdo.FixedReserves(t,j,k);
2124 r.R = z;
2125 r.prc = z;
2126 r.mu = struct('l', z, 'u', z, 'Pmax', z);
2127 r.totalcost = sum(om.get_soln('qdc', 'Rcost', {t,j,k}));
2128 r.R(r.igr) = om.get_soln('var', 'R', {t,j,k}) * baseMVA;
2129 R_mu_l = om.get_soln('lin', 'mu_l', 'Rreq', {t,j,k});
2130 for gg = r.igr
2131 iz = find(r.zones(:, gg));
2132 r.prc(gg) = sum(R_mu_l(iz)) / baseMVA;
2133 end
2134 r.mu.l(r.igr) = om.get_soln('var', 'mu_l', 'R', {t,j,k}) / baseMVA;
2135 r.mu.u(r.igr) = om.get_soln('var', 'mu_u', 'R', {t,j,k}) / baseMVA;
2136 r.mu.Pmax(r.igr) = om.get_soln('lin', 'mu_u', 'Pg_plus_R', {t,j,k}) / baseMVA;
2137 mpc.reserves = r;
2138 end
2139 mdo.flow(t,j,k).mpc = mpc;
2140 end
2141 end
2142
2143 mdo.results.Pc(:,t) = baseMVA * om.get_soln('var', 'Pc', {t});
2144 mdo.results.Rpp(:,t) = baseMVA * om.get_soln('var', 'Rpp', {t});
2145 mdo.results.Rpm(:,t) = baseMVA * om.get_soln('var', 'Rpm', {t});
2146 if ns
2147 mdo.results.Sm(:,t) = baseMVA * om.get_soln('var', 'Sm', {t});
2148 mdo.results.Sp(:,t) = baseMVA * om.get_soln('var', 'Sp', {t});
2149 end
2150 end
2151
2152 for t = 1:mdo.idx.ntramp
2153 mdo.results.Rrp(:,t) = baseMVA * om.get_soln('var', 'Rrp', {t});
2154 mdo.results.Rrm(:,t) = baseMVA * om.get_soln('var', 'Rrm', {t});
2155 end
2156
2157
2158 mdo.results.GenPrices = zeros(ng, nt);
2159 mdo.results.CondGenPrices = zeros(ng, nt);
2160 for t = 1:nt
2161 pp = zeros(ng,1);
2162 for j = 1:mdo.idx.nj(t)
2163 for k = 1:mdo.idx.nc(t,j)+1
2164 pp = pp + mdo.flow(t,j,k).mpc.bus(mdo.flow(t,j,k).mpc.gen(:,GEN_BUS), LAM_P);
2165 end
2166 end
2167 mdo.results.GenPrices(:,t) = pp;
2168 mdo.results.CondGenPrices(:, t) = pp / mdo.StepProb(t);
2169 end
2170
2171 mdo.results.RppPrices = zeros(ng, nt);
2172 mdo.results.RpmPrices = zeros(ng, nt);
2173 for t = 1:nt
2174 mdo.results.RppPrices(:, t) = om.get_soln('var', 'mu_l', 'Rpp', {t}) / baseMVA;
2175 mdo.results.RpmPrices(:, t) = om.get_soln('var', 'mu_l', 'Rpm', {t}) / baseMVA;
2176 for j = 1:mdi.idx.nj(t);
2177 for k = 1:mdi.idx.nc(t,j)+1
2178 ii = find(mdi.flow(t,j,k).mpc.gen(:, GEN_STATUS) > 0);
2179 mdo.results.RppPrices(ii, t) = mdo.results.RppPrices(ii, t) + om.get_soln('lin', 'mu_l', 'dPpRp', {t,j,k}) / baseMVA;
2180 mdo.results.RpmPrices(ii, t) = mdo.results.RpmPrices(ii, t) + om.get_soln('lin', 'mu_l', 'dPmRm', {t,j,k}) / baseMVA;
2181 end
2182 end
2183 end
2184
2185 mdo.results.RrpPrices = zeros(ng, mdo.idx.ntramp);
2186 mdo.results.RrmPrices = zeros(ng, mdo.idx.ntramp);
2187
2188 for t = 1:nt-1
2189 for j1 = 1:mdo.idx.nj(t)
2190 for j2 = 1:mdo.idx.nj(t+1)
2191 if mdi.tstep(t+1).TransMask(j2,j1)
2192 mdo.results.RrpPrices(:, t) = mdo.results.RrpPrices(:, t) + om.get_soln('lin', 'mu_l', 'Rrp', {t,j1,j2}) / baseMVA;
2193 mdo.results.RrmPrices(:, t) = mdo.results.RrmPrices(:, t) + om.get_soln('lin', 'mu_l', 'Rrm', {t,j1,j2}) / baseMVA;
2194 end
2195 end
2196 end
2197 end
2198
2199 if ~mdo.OpenEnded
2200 for j1 = 1:mdo.idx.nj(nt)
2201 mdo.results.RrpPrices(:, nt) = mdo.results.RrpPrices(:, nt) + om.get_soln('lin', 'mu_l', 'Rrp', {nt,j1,1}) / baseMVA;
2202 mdo.results.RrmPrices(:, nt) = mdo.results.RrmPrices(:, nt) + om.get_soln('lin', 'mu_l', 'Rrm', {nt,j1,1}) / baseMVA;
2203 end
2204 end
2205
2206 mdo.results.ExpectedRampCost = zeros(ng, mdo.idx.ntramp+1);
2207
2208 for j = 1:mdi.idx.nj(1)
2209 w = mdo.tstep(1).TransMat(j,1);
2210 mdo.results.ExpectedRampCost(:, 1) = mdo.results.ExpectedRampCost(:, 1) ...
2211 + 0.5 * w * mdo.RampWearCostCoeff(:,1) .* (mdo.flow(1,j,1).mpc.gen(:,PG) - mdo.InitialPg).^2;
2212 end
2213
2214 for t = 2:nt
2215 for j2 = 1:mdo.idx.nj(t)
2216 for j1 = 1:mdo.idx.nj(t-1)
2217 w = mdo.tstep(t).TransMat(j2,j1) * mdo.CostWeights(1, j1, t-1);
2218 mdo.results.ExpectedRampCost(:, t) = mdo.results.ExpectedRampCost(:, t) ...
2219 + 0.5 * w * mdo.RampWearCostCoeff(:,t) .* (mdo.flow(t,j2,1).mpc.gen(:,PG) - mdo.flow(t-1,j1,1).mpc.gen(:,PG)) .^2;
2220 end
2221 end
2222 end
2223
2224 if ~mdo.OpenEnded
2225 for j = 1:mdi.idx.nj(nt)
2226 w = mdi.tstep(t+1).TransMat(1, j) * mdi.CostWeights(1, j, nt);
2227 mdo.results.ExpectedRampCost(:, nt+1) = 0.5 * w * mdo.RampWearCostCoeff(:,nt+1) .* (mdo.TerminalPg - mdo.flow(nt,j,1).mpc.gen(:,PG)) .^2;
2228 end
2229 end
2230
2231
2232 mdo.results.ExpectedDispatch = zeros(ng, nt);
2233 for t = 1:nt
2234 pp = sum(mdo.CostWeights(1,1:mdo.idx.nj(t),t)');
2235 for j = 1:mdo.idx.nj(t)
2236 mdo.results.ExpectedDispatch(:,t) = mdo.results.ExpectedDispatch(:,t) + ...
2237 mdo.CostWeights(1,j,t)/pp * mdo.flow(t,j,1).mpc.gen(:,PG);
2238 end
2239 end
2240
2241 if ns && mdo.Storage.ForceCyclicStorage
2242 mdo.Storage.InitialStorage = baseMVA * om.get_soln('var', 'S0');
2243 end
2244
2245 mdo.Storage.ExpectedStorageState = zeros(ns,nt);
2246 if ns
2247 for t = 1:nt
2248 pp = sum(mdo.CostWeights(1,1:mdo.idx.nj(t),t)');
2249 for j = 1:mdo.idx.nj(t)
2250 Lfj = mdo.tstep(t).Lf( j:mdo.idx.nj(t):(ns-1)*mdo.idx.nj(t)+j, :);
2251 Ngj = mdo.tstep(t).Ng( j:mdo.idx.nj(t):(ns-1)*mdo.idx.nj(t)+j, :);
2252 Nhj = mdo.tstep(t).Nh( j:mdo.idx.nj(t):(ns-1)*mdo.idx.nj(t)+j, :);
2253 mdo.Storage.ExpectedStorageState(:,t) = ...
2254 mdo.Storage.ExpectedStorageState(:,t) + ...
2255 baseMVA * mdo.CostWeights(1,j,t)/pp * ...
2256 ( Lfj * mdo.Storage.InitialStorage/baseMVA + ...
2257 (Ngj + Nhj) * mdo.QP.x );
2258 end
2259 end
2260 mdo.Storage.ExpectedStorageDispatch = ...
2261 mdo.results.ExpectedDispatch(mdo.Storage.UnitIdx, :);
2262 end
2263
2264
2265 if ntds
2266 if nzds
2267 mdo.results.Z = zeros(nzds, ntds);
2268 for t = 1:ntds
2269 mdo.results.Z(:,t) = om.get_soln('var', 'Z', {t});
2270 end
2271 end
2272 mdo.results.Y = zeros(nyds, ntds);
2273 if nyds
2274 for t = 1:ntds
2275 AA = om.params_lin_constraint('DSy', {t});
2276 mdo.results.Y(:, t) = AA * mdo.QP.x;
2277 end
2278 end
2279 end
2280 mdo.results.f = mdo.QP.f;
2281 end
2282 mdo.results.success = success;
2283 mdo.results.SolveTime = toc(t0);
2284 end
2285
2286 mdo.results.SetupTime = et_setup;
2287
2288 if verbose
2289 fprintf('- MOST: Done.\n\n');
2290 end