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