0001 function t_most_w_ds(quiet, solver, verbose)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 if nargin < 3
0014 verbose = 0;
0015 if nargin < 2
0016 solver = '';
0017 if nargin < 1
0018 quiet = 0;
0019 end
0020 end
0021 end
0022
0023 include_MIPS = 0;
0024
0025 n_tests = 2;
0026
0027 t_begin(n_tests, quiet);
0028
0029 casefile = 'c118swf';
0030 solnfile = 't_most_w_ds_z';
0031
0032
0033 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0034 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0035 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0036 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0037 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0038 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0039 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0040 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0041 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0042 [CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
0043 CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
0044 CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
0045 CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
0046 CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
0047 CT_MODCOST_X] = idx_ct;
0048
0049 if have_feature('cplex') || have_feature('gurobi') || ...
0050 have_feature('mosek') || have_feature('quadprog_ls') || include_MIPS
0051 mdi = md_init;
0052
0053
0054 if isempty(solver) || strcmp(upper(solver), 'DEFAULT')
0055 if have_feature('mosek')
0056 solver = 'MOSEK';
0057 elseif have_feature('cplex')
0058 solver = 'CPLEX';
0059 elseif have_feature('gurobi')
0060 solver = 'GUROBI';
0061 elseif have_feature('quadprog_ls')
0062 solver = 'OT';
0063 else
0064 solver = 'MIPS';
0065 end
0066 end
0067 mpopt = mpoption('most.solver', solver, 'verbose', verbose);
0068
0069
0070 if have_feature('cplex')
0071 mpopt = mpoption(mpopt, 'cplex.opts.threads', 2);
0072 end
0073 if have_feature('gurobi')
0074 mpopt = mpoption(mpopt, 'gurobi.method', 2);
0075 mpopt = mpoption(mpopt, 'gurobi.threads', 2);
0076 mpopt = mpoption(mpopt, 'gurobi.opts.BarConvTol', 1e-6);
0077 mpopt = mpoption(mpopt, 'gurobi.opts.FeasibilityTol', 1e-4);
0078 mpopt = mpoption(mpopt, 'gurobi.opts.OptimalityTol', 1e-5);
0079 end
0080 if have_feature('quadprog_ls')
0081 mpopt = mpoption(mpopt, 'quadprog.TolFun', 1e-13);
0082 end
0083 if have_feature('mosek')
0084 mpopt = mpoption(mpopt, 'mosek.num_threads', 2);
0085 else
0086 mpopt = mpoption(mpopt, 'mips.max_it', 500);
0087 if have_feature('pardiso')
0088 mpopt = mpoption(mpopt, 'mips.linsolver', 'PARDISO');
0089 end
0090 end
0091
0092
0093
0094
0095
0096
0097
0098
0099 mdi.mpc = loadcase(casefile);
0100 mdi.InitialPg = mdi.mpc.gen(:,PG);
0101 nt = 24;
0102 ng = size(mdi.mpc.gen, 1);
0103 mdi.idx.nt = nt;
0104 PositiveActiveReservePrice = ones(ng,1);
0105 PositiveActiveReserveQuantity = 0.25*mdi.mpc.gen(:,PMAX);
0106 NegativeActiveReservePrice = ones(ng,1);
0107 NegativeActiveReserveQuantity = PositiveActiveReserveQuantity;
0108 PositiveActiveDeltaPrice = ones(ng,1);
0109 NegativeActiveDeltaPrice = ones(ng,1);
0110 PositiveLoadFollowReservePrice = ones(ng,1);
0111 PositiveLoadFollowReserveQuantity = 0.5*mdi.mpc.gen(:,PMAX);
0112 NegativeLoadFollowReservePrice = ones(ng,1);
0113 NegativeLoadFollowReserveQuantity = PositiveLoadFollowReserveQuantity;
0114
0115
0116
0117 mdi.mpc.gen(:,RAMP_10) = 1.0 * mdi.mpc.gen(:,PMAX);
0118 mdi.mpc.gen(:,RAMP_AGC) = 1.0 * mdi.mpc.gen(:,PMAX);
0119 mdi.mpc.gen(:,RAMP_30) = 1.0 * mdi.mpc.gen(:,PMAX);
0120
0121 mdi.RampWearCostCoeff = 0.05 * ones(ng,1);
0122 for t = 2:nt
0123 mdi.RampWearCostCoeff(:, t) = mdi.RampWearCostCoeff(:, 1);
0124 end
0125 mdi.Storage(1).UnitIdx = mdi.mpc.iess;
0126 ns = length(mdi.Storage.UnitIdx);
0127 Minstor = zeros(ns,1);
0128 Maxstor = 200 * ones(ns,1);
0129
0130
0131 mdi.Storage.InitialStorage = 50 * ones(ns,1);
0132 mdi.Storage.InitialStorageLowerBound = 50*ones(ns,1);
0133 mdi.Storage.InitialStorageUpperBound = 50*ones(ns,1);
0134 mdi.Storage.OutEff = 0.95 * ones(ns,1);
0135 mdi.Storage.InEff = 0.9 * ones(ns ,1);
0136 mdi.Storage.InitialStorageCost = 35 * ones(ns, 1);
0137 mdi.Storage.TerminalStoragePrice = 35 * ones(ns, 1);
0138 mdi.Storage.TerminalChargingPrice0 = 35 * ones(ns, 1);
0139 mdi.Storage.TerminalDischargingPrice0 = 35 * ones(ns, 1);
0140 mdi.Storage.TerminalChargingPriceK = 10 * ones(ns, 1);
0141 mdi.Storage.TerminalDischargingPriceK = 40 * ones(ns, 1);
0142 mpopt = mpoption(mpopt, 'most.storage.terminal_target', 0);
0143 mdi.Storage.ExpectedTerminalStorageAim = mdi.Storage.InitialStorage;
0144 mdi.Storage.LossFactor = zeros(ns,1);
0145 mdi.Storage.IncludeValueOfTerminalStorage = 1;
0146 mpopt = mpoption(mpopt, 'most.storage.cyclic', 1);
0147
0148 for t = 1:nt
0149 mdi.offer(t).gencost = mdi.mpc.gencost;
0150 mdi.offer(t).PositiveActiveReservePrice = PositiveActiveReservePrice;
0151 mdi.offer(t).PositiveActiveReserveQuantity = PositiveActiveReserveQuantity;
0152 mdi.offer(t).NegativeActiveReservePrice = NegativeActiveReservePrice;
0153 mdi.offer(t).NegativeActiveReserveQuantity = NegativeActiveReserveQuantity;
0154 mdi.offer(t).PositiveActiveDeltaPrice = PositiveActiveDeltaPrice;
0155 mdi.offer(t).NegativeActiveDeltaPrice = NegativeActiveDeltaPrice;
0156 mdi.offer(t).PositiveLoadFollowReservePrice = PositiveLoadFollowReservePrice;
0157 mdi.offer(t).PositiveLoadFollowReserveQuantity = PositiveLoadFollowReserveQuantity;
0158 mdi.offer(t).NegativeLoadFollowReservePrice = NegativeLoadFollowReservePrice;
0159 mdi.offer(t).NegativeLoadFollowReserveQuantity = NegativeLoadFollowReserveQuantity;
0160 mdi.Storage.MinStorageLevel(:,t) = Minstor;
0161 mdi.Storage.MaxStorageLevel(:,t) = Maxstor;
0162 end
0163 mdi.Storage.MinStorageLevel(:,nt+1) = Minstor;
0164 mdi.Storage.MaxStorageLevel(:,nt+1) = Maxstor;
0165
0166
0167
0168
0169 mdi.UC.CommitSched = ones(ng,nt);
0170
0171 mdi.Delta_T = 1;
0172
0173 loadprof = [ 0.6 0.6 0.6 0.7 0.75 0.8 0.9 1.1 1.2 1.3 1.4 1.4 1.3 1.4 1.4 1.4 1.4 1.3 1.1 1.1 1.0 0.9 0.8 0.7 ];
0174
0175
0176
0177 partialcontabrow =[ 1 0 CT_TBUS 0 PD CT_REL ];
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187 for t = 1:nt
0188 mdi.tstep(t).OpCondSched(1).tab = [ partialcontabrow loadprof(t) ];
0189 end
0190
0191
0192 for t = 1:nt
0193 mdi.tstep(t).OpCondSched(2).tab = mdi.tstep(t).OpCondSched(1).tab;
0194 mdi.tstep(t).OpCondSched(2).tab(1,7) = 1.1*mdi.tstep(t).OpCondSched(2).tab(1,7);
0195 mdi.tstep(t).OpCondSched(3).tab = mdi.tstep(t).OpCondSched(1).tab;
0196 mdi.tstep(t).OpCondSched(3).tab(1,7) = 0.9*mdi.tstep(t).OpCondSched(1).tab(1,7);
0197 end
0198
0199
0200 contab = [
0201 1 0.01 CT_TGEN 2 GEN_STATUS CT_REP 0 ;
0202 2 0.01 CT_TGEN 5 GEN_STATUS CT_REP 0 ;
0203 ];
0204
0205 for t = 1:nt
0206 for j = 1:3
0207 mdi.cont(t,j).contab = contab;
0208 end
0209 end
0210
0211
0212
0213 mdi.tstep(1).TransMat = [ 1/3;
0214 1/3
0215 1/3];
0216 for t = 2:nt
0217 mdi.tstep(t).TransMat = 1/3 * ones(3,3);
0218 end
0219
0220
0221 ntds = 24;
0222 mdi.idx.ntds = ntds;
0223 m1 = 8;
0224 m2 = 12;
0225 B = sparse(m1*m2, ng);
0226 ilist = [ 2 3 4 5 3 4 5 6 3 4 5 7 3 4 5 7 3 4 5 6 4 5 6 7 2 3 4 ];
0227 jlist = [ 2 2 2 2 3 3 3 3 5 5 5 5 6 6 6 6 7 7 7 7 8 8 8 8 11 11 11 ];
0228 for i = 1:length(mdi.mpc.icoal)
0229 B((jlist(i)-1)*m1+ilist(i), mdi.mpc.icoal(i)) = 0.1;
0230 end
0231 A = mkdif(m1, m2, 0.5, 0.97, [1.0 0]);
0232 C = [];
0233 D = [];
0234 zmin = zeros(m1*m2, 1);
0235 zmax = 100*ones(m1*m2, 1);
0236 ymin = 0;
0237 ymax = 100;
0238 for t = 1:ntds
0239 mdi.dstep(t).A = A;
0240 mdi.dstep(t).B = B;
0241 mdi.dstep(t).C = C;
0242 mdi.dstep(t).D = D;
0243 mdi.dstep(t).zmin = zmin;
0244 mdi.dstep(t).zmax = zmax;
0245 mdi.dstep(t).ymin = ymin;
0246 mdi.dstep(t).ymax = ymax;
0247 end
0248 mdi.z1 = zeros(m1*m2, 1);
0249
0250 mdo = most(mdi, mpopt);
0251
0252 s = load(solnfile);
0253
0254 t = 'objective function value (f)';
0255 t_is(mdo.QP.f, 1575531.9, -0.5, t);
0256
0257
0258
0259
0260
0261 t = 'dynamical system state (Z)';
0262 t_is(mdo.results.Z, s.Z, 3.7, t);
0263 else
0264 t_skip(2, 'requires MOSEK, CPLEX, Gurobi or quadprog');
0265 end
0266
0267
0268
0269
0270
0271
0272 t_end
0273
0274 function A = mkdif(m1, m2, alpha, r, w)
0275
0276
0277
0278
0279
0280
0281
0282
0283 if norm(w) > 1
0284 w = w / norm(w);
0285 end
0286
0287 n = m1 * m2;
0288 A = sparse(n, n);
0289 north = alpha / 4;
0290 east = north;
0291 south = north;
0292 west = north;
0293 self = (1-alpha);
0294 if w(1) > 0
0295 west = west + w(1)*self;
0296 self = (1-w(1)) * self;
0297 elseif w(2) < 0
0298 east = east + (-w(1))*self;
0299 self = (1+w(1)) * self;
0300 end
0301 if w(2) > 0
0302 south = south + w(2)*self;
0303 self = (1-w(2)) * self;
0304 elseif w(2) < 0
0305 north = north + (-w(2))*self;
0306 self = (1+w(2)) * self;
0307 end
0308 tot = self + north + south + east + west;
0309 self = (r/tot) * self;
0310 north = (r/tot) * north;
0311 south = (r/tot) * south;
0312 east = (r/tot) * east;
0313 west = (r/tot) * west;
0314
0315 for i = 1:m1
0316 for j = 1:m2
0317 A((j-1)*m1+i, (j-1)*m1+i) = self;
0318
0319 if i > 1
0320 A((j-1)*m1+i, (j-1)*m1+i-1) = north;
0321 end
0322
0323 if i < m1
0324 A((j-1)*m1+i, (j-1)*m1+i+1) = south;
0325 end
0326
0327 if j > 1
0328 A((j-1)*m1+i, (j-2)*m1+i) = west;
0329 end
0330
0331 if j < m2
0332 A((j-1)*m1+i, j*m1+i) = east;
0333 end
0334 end
0335 end