Description of t_most_w

0001 function t_most_w_ds(quiet, solver, verbose)
0002 %T_MOST_W_DS  Test for MOST with dynamical system constraints.
0003 
0004 %   MOST
0005 %   Copyright (c) 2015-2016, Power Systems Engineering Research Center (PSERC)
0006 %   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
0007 %   and Ray Zimmerman, PSERC Cornell
0008 %
0009 %   This file is part of MOST.
0010 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0011 %   See https://github.com/MATPOWER/most for more info.
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;   %% set to 1, to attempt even if MIPS is the best solver
0024                     %% available (takes a LONG time and currently fails)
0025 n_tests = 1;
0026 
0027 t_begin(n_tests, quiet);
0028 
0029 casefile = 'c118swf';
0030 solnfile = 't_most_w_ds_z';
0031 
0032 %% define named indices into data matrices
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_fcn('cplex') || have_fcn('gurobi') || have_fcn('mosek') || ...
0050         have_fcn('quadprog_ls') || include_MIPS
0051     mdi = md_init;
0052 
0053     %% choose solver
0054     if isempty(solver) || strcmp(upper(solver), 'DEFAULT')
0055         if have_fcn('mosek')
0056             solver = 'MOSEK';
0057         elseif have_fcn('cplex')
0058             solver = 'CPLEX';
0059         elseif have_fcn('gurobi')
0060             solver = 'GUROBI';
0061         elseif have_fcn('quadprog_ls')
0062             solver = 'OT';
0063         else
0064             solver = 'MIPS';
0065         end
0066     end
0067     mpopt = mpoption('most.solver', solver, 'verbose', verbose);
0068 
0069     %% set options
0070     if have_fcn('cplex')
0071         mpopt = mpoption(mpopt, 'cplex.opts.threads', 2);   % set this manually here
0072     end
0073     if have_fcn('gurobi')
0074         mpopt = mpoption(mpopt, 'gurobi.method', 2);        %% barrier
0075         mpopt = mpoption(mpopt, 'gurobi.threads', 2);
0076         mpopt = mpoption(mpopt, 'gurobi.opts.BarConvTol', 1e-6);        %% 1e-8
0077         mpopt = mpoption(mpopt, 'gurobi.opts.FeasibilityTol', 1e-4);    %% 1e-6
0078         mpopt = mpoption(mpopt, 'gurobi.opts.OptimalityTol', 1e-5);     %% 1e-6
0079     end
0080     if have_fcn('quadprog_ls')
0081         mpopt = mpoption(mpopt, 'quadprog.TolFun', 1e-13);
0082     end
0083     if have_fcn('mosek')
0084         mpopt = mpoption(mpopt, 'mosek.num_threads', 2);
0085     else
0086         mpopt = mpoption(mpopt, 'mips.max_it', 500);
0087         if have_fcn('pardiso')
0088             mpopt = mpoption(mpopt, 'mips.linsolver', 'PARDISO');
0089         end
0090     end
0091     %% use e.g. t_most_w_ds(0, 'MOSEK') instead of uncommenting these lines
0092     % mpopt = mpoption(mpopt, 'most.solver', 'MOSEK');
0093     % mpopt = mpoption(mpopt, 'most.solver', 'CPLEX');
0094     % mpopt = mpoption(mpopt, 'most.solver', 'GUROBI');
0095     % mpopt = mpoption(mpopt, 'most.solver', 'OT');
0096     % mpopt = mpoption(mpopt, 'most.solver', 'BPMPD');
0097     % mpopt = mpoption(mpopt, 'most.solver', 'MIPS');
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     %mdi.mpc.gen(:,RAMP_10) = 0.20 * mdi.mpc.gen(PMAX);
0115     %mdi.mpc.gen(:,RAMP_AGC) = 0.20 * mdi.mpc.gen(PMAX);
0116     %mdi.mpc.gen(:,RAMP_30) = 0.50 * mdi.mpc.gen(PMAX);
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);   % (i, t) note different scheme!
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     %mdi.Storage.MinStorageLevel    = zeros(ns,1);
0130     %mdi.Storage.MaxStorageLevel    = 200 * ones(ns,1);
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); % applied to psc_tij0, psd_tij0 (non-terminal states)
0138     mdi.Storage.TerminalChargingPrice0     = 35 * ones(ns, 1); % applied to psc_tijk (contingency terminal states)
0139     mdi.Storage.TerminalDischargingPrice0  = 35 * ones(ns, 1); % applied to psd_tijk (contingency terminal states)
0140     mdi.Storage.TerminalChargingPriceK     = 10 * ones(ns, 1); % applied to psc_tij0 (end-of-horizon terminal states)
0141     mdi.Storage.TerminalDischargingPriceK  = 40 * ones(ns, 1); % applied to psd_tij0 (end-of-horizon terminal states)
0142     mpopt = mpoption(mpopt, 'most.storage.terminal_target', 0);
0143     mdi.Storage.ExpectedTerminalStorageAim = mdi.Storage.InitialStorage;  % expected terminal storage if mpopt.most.storage.terminal_target is true
0144     mdi.Storage.LossFactor         = zeros(ns,1);  % fraction of storage lost in each period
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;  % Needed if mpopt.most.storage.cyclic
0164     mdi.Storage.MaxStorageLevel(:,nt+1) = Maxstor;
0165 
0166     %mdi.Storage.MinStorageLevel(:,4) = [10 ; 10];  % is this enough to create infeasibility?
0167     %mdi.Storage.MaxStorageLevel(:,4) = [10; 50 ];
0168 
0169     mdi.UC.CommitSched = ones(ng,nt);
0170 
0171     mdi.Delta_T = 1;
0172     %            0:0  1:00 2:00 3:00 4:00  5:00  6:00  7:00 8:00 9:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00 23:00
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     %                label  probability   type      row     column      chg type  newvalue
0177     partialcontabrow =[ 1       0        CT_TBUS     0        PD         CT_REL ];
0178     %mdi.tstep(1).OpCondSched(1).tab= [
0179     %                   1        0        CT_TBUS     0        PD         CT_REL     0.8  ];
0180     %mdi.tstep(2).OpCondSched(1).tab= [
0181     %                   1        0        CT_TBUS     0        PD         CT_REL     1.0  ];
0182     %mdi.tstep(3).OpCondSched(1).tab= [
0183     %                   1        0        CT_TBUS     0        PD         CT_REL     1.2  ];
0184     %mdi.tstep(4).OpCondSched(1).tab= [
0185     %                   1        0        CT_TBUS     0        PD         CT_REL     0.9 ];
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 = [%         1       0.01      CT_TBUS     0        PD         CT_REL     1.05 ;
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  % mdi.idx.nj(t)
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 = 'dynamical system state (Z)';
0255     t_is(mdo.results.Z, s.Z, 3.7, t);
0256 else
0257     t_skip(1, 'requires MOSEK, CPLEX, Gurobi or quadprog');
0258 end
0259 
0260 % YorN = input('Play movie? (y/n) : ', 's');
0261 % if strcmp(upper(YorN(1)), 'Y')
0262 %     domovie;
0263 % end
0264 
0265 t_end
0266 
0267 function A = mkdif(m1, m2, alpha, r, w)
0268 % A = mkdif(m1, m2, r, w)
0269 %
0270 % computes an A matrix for the difussion equations in an m1 x m2 grid,
0271 % using a difussion speed factor alpha <= 1, a dissipation factor r < 1 and
0272 % a "wind" vector w = [wx wy] that roughly tells where the pollutants go.
0273 
0274 % Carlos Murillo.  As naive as can be.
0275 
0276 if norm(w) > 1
0277   w = w / norm(w);
0278 end
0279 
0280 n = m1 * m2;
0281 A = sparse(n, n);
0282 north = alpha / 4;
0283 east = north;
0284 south = north;
0285 west = north;
0286 self = (1-alpha);
0287 if w(1) > 0
0288   west = west + w(1)*self;
0289   self = (1-w(1)) * self;
0290 elseif w(2) < 0
0291   east = east + (-w(1))*self;
0292   self = (1+w(1)) * self;
0293 end
0294 if w(2) > 0
0295   south = south + w(2)*self;
0296   self = (1-w(2)) * self;
0297 elseif w(2) < 0
0298   north = north + (-w(2))*self;
0299   self = (1+w(2)) * self;
0300 end
0301 tot = self + north + south + east + west;
0302 self = (r/tot) * self;
0303 north = (r/tot) * north;
0304 south = (r/tot) * south;
0305 east = (r/tot) * east;
0306 west = (r/tot) * west;
0307 
0308 for i = 1:m1
0309   for j = 1:m2
0310     A((j-1)*m1+i, (j-1)*m1+i) = self;
0311     % North
0312     if i > 1
0313       A((j-1)*m1+i, (j-1)*m1+i-1) = north;
0314     end
0315     % South
0316     if i < m1
0317       A((j-1)*m1+i, (j-1)*m1+i+1) = south;
0318     end
0319     % West
0320     if j > 1
0321       A((j-1)*m1+i, (j-2)*m1+i) = west;
0322     end
0323     % East
0324     if j < m2
0325       A((j-1)*m1+i, j*m1+i) = east;
0326     end
0327   end
0328 end
t_most_w_ds

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
