Description of t_most_w

0001 function t_most_w_ds(quiet)
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 http://www.pserc.cornell.edu/matpower/ for more info.
0012 
0013 if nargin < 1
0014     quiet = 0;
0015 end
0016 
0017 include_MIPS = 0;   %% set to 1, to attempt even if MIPS is the best solver
0018                     %% available (takes a LONG time and currently fails)
0019 n_tests = 1;
0020 
0021 t_begin(n_tests, quiet);
0022 
0023 casefile = 'c118swf';
0024 solnfile = 't_most_w_ds_z';
0025 
0026 %% define named indices into data matrices
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);   % set this manually here
0052     elseif have_fcn('gurobi')
0053         mpopt = mpoption(mpopt, 'most.solver', 'GUROBI');
0054         mpopt = mpoption(mpopt, 'gurobi.method', 2);        %% barrier
0055         mpopt = mpoption(mpopt, 'gurobi.threads', 2);
0056         mpopt = mpoption(mpopt, 'gurobi.opts.BarConvTol', 1e-7);        %% 1e-8
0057         mpopt = mpoption(mpopt, 'gurobi.opts.FeasibilityTol', 1e-5);    %% 1e-6
0058         mpopt = mpoption(mpopt, 'gurobi.opts.OptimalityTol', 1e-5);     %% 1e-6
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     % mpopt = mpoption(mpopt, 'most.solver', 'CLP');
0073     % mpopt = mpoption(mpopt, 'most.solver', 'IPOPT');
0074     % mpopt = mpoption(mpopt, 'most.solver', 'MIPS');
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     %mdi.mpc.gen(:,RAMP_10) = 0.20 * mdi.mpc.gen(PMAX);
0092     %mdi.mpc.gen(:,RAMP_AGC) = 0.20 * mdi.mpc.gen(PMAX);
0093     %mdi.mpc.gen(:,RAMP_30) = 0.50 * mdi.mpc.gen(PMAX);
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);   % (i, t) note different scheme!
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     %mdi.Storage.MinStorageLevel    = zeros(ns,1);
0107     %mdi.Storage.MaxStorageLevel    = 200 * ones(ns,1);
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); % applied to psc_tij0, psd_tij0 (non-terminal states)
0115     mdi.Storage.TerminalChargingPrice0     = 35 * ones(ns, 1); % applied to psc_tijk (contingency terminal states)
0116     mdi.Storage.TerminalDischargingPrice0  = 35 * ones(ns, 1); % applied to psd_tijk (contingency terminal states)
0117     mdi.Storage.TerminalChargingPriceK     = 10 * ones(ns, 1); % applied to psc_tij0 (end-of-horizon terminal states)
0118     mdi.Storage.TerminalDischargingPriceK  = 40 * ones(ns, 1); % applied to psd_tij0 (end-of-horizon terminal states)
0119     mpopt = mpoption(mpopt, 'most.storage.terminal_target', 0);
0120     mdi.Storage.ExpectedTerminalStorageAim = mdi.Storage.InitialStorage;  % expected terminal storage if mpopt.most.storage.terminal_target is true
0121     mdi.Storage.LossFactor         = zeros(ns,1);  % fraction of storage lost in each period
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;  % Needed if mpopt.most.storage.cyclic
0141     mdi.Storage.MaxStorageLevel(:,nt+1) = Maxstor;
0142 
0143     %mdi.Storage.MinStorageLevel(:,4) = [10 ; 10];  % is this enough to create infeasibility?
0144     %mdi.Storage.MaxStorageLevel(:,4) = [10; 50 ];
0145 
0146     mdi.UC.CommitSched = ones(ng,nt);
0147 
0148     mdi.Delta_T = 1;
0149     %            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
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     %                label  probability   type      row     column      chg type  newvalue
0154     partialcontabrow =[ 1       0        CT_TBUS     0        PD         CT_REL ];
0155     %mdi.tstep(1).OpCondSched(1).tab= [
0156     %                   1        0        CT_TBUS     0        PD         CT_REL     0.8  ];
0157     %mdi.tstep(2).OpCondSched(1).tab= [
0158     %                   1        0        CT_TBUS     0        PD         CT_REL     1.0  ];
0159     %mdi.tstep(3).OpCondSched(1).tab= [
0160     %                   1        0        CT_TBUS     0        PD         CT_REL     1.2  ];
0161     %mdi.tstep(4).OpCondSched(1).tab= [
0162     %                   1        0        CT_TBUS     0        PD         CT_REL     0.9 ];
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 = [%         1       0.01      CT_TBUS     0        PD         CT_REL     1.05 ;
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  % mdi.idx.nj(t)
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 % YorN = input('Play movie? (y/n) : ', 's');
0238 % if strcmp(upper(YorN(1)), 'Y')
0239 %     domovie;
0240 % end
0241 
0242 t_end
0243 
0244 function A = mkdif(m1, m2, alpha, r, w)
0245 % A = mkdif(m1, m2, r, w)
0246 %
0247 % computes an A matrix for the difussion equations in an m1 x m2 grid,
0248 % using a difussion speed factor alpha <= 1, a dissipation factor r < 1 and
0249 % a "wind" vector w = [wx wy] that roughly tells where the pollutants go.
0250 
0251 % Carlos Murillo.  As naive as can be.
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     % North
0289     if i > 1
0290       A((j-1)*m1+i, (j-1)*m1+i-1) = north;
0291     end
0292     % South
0293     if i < m1
0294       A((j-1)*m1+i, (j-1)*m1+i+1) = south;
0295     end
0296     % West
0297     if j > 1
0298       A((j-1)*m1+i, (j-2)*m1+i) = west;
0299     end
0300     % East
0301     if j < m2
0302       A((j-1)*m1+i, j*m1+i) = east;
0303     end
0304   end
0305 end
t_most_w_ds

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
