plot_storage

PLOT_STORAGE   Plot storage unit results

   PLOT_STORAGE(MD)
   PLOT_STORAGE(MD, IDX)
   PLOT_STORAGE(MD, IDX, '<option1_name>', '<option1_value', ...)
   PLOT_STORAGE(MD, IDX, OPT)
   PLOT_STORAGE(MD, IDX, OPT, '<option1_name>', '<option1_value', ...)
   H = PLOT_STORAGE(MD, ...)

   IDX is the storage unit index. If IDX is a vector, it sums them first,
   if empty, it includes all storage units in MD. Options can be
   specified as an OPT struct or as individual pairs of 'name' and 'value'
   arguments. The possible options include the following, where the default
   is shown in parenthesis:
       'saveit'        (false) flag to indicate whether to create PDF file
       'saveall'       (false) flag to indicate whether to create individual
                       PDF files for each element of IDX, as well as
                       aggregate, when IDX is a vector (or empty)
       'savepath'      ('') path to directory to save files in
       'savename'      ('stored-energy-%s.pdf') name of PDF file
                       %s is optional placeholder for storage unit index
       'separation'    (0.8)  %% separation of beginning/end of period (0-1)
           0   = beginning & end of period t, staircase (dispatches are
                 vertical lines at t)
           0.5 = evenly separated (both dispatch and transitions visible)
           1   = beginning of t aligned with end of t-1, smooth,
                 (transitions are vertical displacements at t+/-0.5)
       'sort_tol'      (1e-6) round to nearest sort_tol for sorting
       'size_factor'   (1) to scale font/marker sizes in case you want to
                       do sub-plots or something
       'show_grid'                 (true) vertical lines to divide periods
       'show_expected_initial'     (true)
       'show_expected_final'       (true)
       'show_adjusted_dispatches'  (true)
       'show_dispatches'           (false)

   Returns handle to current figure window.

0001 function hh = plot_storage(md, idx, varargin)
0002 %PLOT_STORAGE   Plot storage unit results
0003 %
0004 %   PLOT_STORAGE(MD)
0005 %   PLOT_STORAGE(MD, IDX)
0006 %   PLOT_STORAGE(MD, IDX, '<option1_name>', '<option1_value', ...)
0007 %   PLOT_STORAGE(MD, IDX, OPT)
0008 %   PLOT_STORAGE(MD, IDX, OPT, '<option1_name>', '<option1_value', ...)
0009 %   H = PLOT_STORAGE(MD, ...)
0010 %
0011 %   IDX is the storage unit index. If IDX is a vector, it sums them first,
0012 %   if empty, it includes all storage units in MD. Options can be
0013 %   specified as an OPT struct or as individual pairs of 'name' and 'value'
0014 %   arguments. The possible options include the following, where the default
0015 %   is shown in parenthesis:
0016 %       'saveit'        (false) flag to indicate whether to create PDF file
0017 %       'saveall'       (false) flag to indicate whether to create individual
0018 %                       PDF files for each element of IDX, as well as
0019 %                       aggregate, when IDX is a vector (or empty)
0020 %       'savepath'      ('') path to directory to save files in
0021 %       'savename'      ('stored-energy-%s.pdf') name of PDF file
0022 %                       %s is optional placeholder for storage unit index
0023 %       'separation'    (0.8)  %% separation of beginning/end of period (0-1)
0024 %           0   = beginning & end of period t, staircase (dispatches are
0025 %                 vertical lines at t)
0026 %           0.5 = evenly separated (both dispatch and transitions visible)
0027 %           1   = beginning of t aligned with end of t-1, smooth,
0028 %                 (transitions are vertical displacements at t+/-0.5)
0029 %       'sort_tol'      (1e-6) round to nearest sort_tol for sorting
0030 %       'size_factor'   (1) to scale font/marker sizes in case you want to
0031 %                       do sub-plots or something
0032 %       'show_grid'                 (true) vertical lines to divide periods
0033 %       'show_expected_initial'     (true)
0034 %       'show_expected_final'       (true)
0035 %       'show_adjusted_dispatches'  (true)
0036 %       'show_dispatches'           (false)
0037 %
0038 %   Returns handle to current figure window.
0039 
0040 %   TO DO: Do initial invisible plot to get v axis parameters.
0041 
0042 %   MOST
0043 %   Copyright (c) 2013-2020, Power Systems Engineering Research Center (PSERC)
0044 %   by Ray Zimmerman, PSERC Cornell
0045 %
0046 %   This file is part of MOST.
0047 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0048 %   See https://github.com/MATPOWER/most for more info.
0049 
0050 %% define named indices into data matrices
0051 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0052     MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0053     QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0054 
0055 %% initialize some variables
0056 my_xlabel = 'Period';
0057 nt = md.idx.nt;
0058 ns = md.idx.ns;
0059 nj_max = max(md.idx.nj);
0060 baseMVA = md.mpc.baseMVA;
0061 
0062 %% input args
0063 if nargin < 2
0064     idx = [];
0065 end
0066 if isempty(idx)
0067     idx = (1:ns)';
0068 end
0069 nidx = length(idx);
0070 if nidx > 1 && size(idx, 1) == 1
0071     idx = idx';     %% convert row vector to column vector
0072 end
0073 g = md.Storage.UnitIdx(idx);
0074 b = md.mpc.gen(g, GEN_BUS);
0075 
0076 %% default options
0077 opt = struct( ...
0078     'saveit', false, ...
0079     'saveall', false, ...
0080     'savepath', '', ...
0081     'savename', 'stored-energy-%s.pdf', ...
0082     'separation', 0.8, ...
0083     'sort_tol', 1e-6, ...
0084     'size_factor', 1, ...
0085     'show_grid', true, ...
0086     'show_adjusted_dispatches', true, ...
0087     'show_expected_initial', true, ...
0088     'show_expected_final', true, ...
0089     'show_dispatches', false );
0090 
0091 %% process options
0092 if mod(length(varargin), 2) %% odd number of options, first must be OPT struct
0093     if ~isstruct(varargin{1})
0094         error('plot_storage: Single OPT argument must be a struct');
0095     end
0096     myopt = varargin{1};
0097     k = 2;
0098 else                        %% even number of options
0099     myopt = struct;
0100     k = 1;
0101 end
0102 while k < length(varargin)
0103     opt_name = varargin{k};
0104     opt_val  = varargin{k+1};
0105     if ~isfield(opt, opt_name)
0106         error('plot_name: ''%s'' is not a valid option name', opt_name);
0107     end
0108     myopt.(opt_name) = opt_val;
0109     k = k + 2;
0110 end
0111 fields = fieldnames(myopt);
0112 for f = 1:length(fields)
0113     opt.(fields{f}) = myopt.(fields{f});
0114 end
0115 
0116 %% call recursively for individual plots if indicated
0117 if opt.saveall && nidx > 1
0118     if isempty(strfind(opt.savename, '%s'))
0119         error('plot_storage: ''savename'' must include a ''%%s'' placeholder when ''saveall'' option is true.');
0120     end
0121     for i = 1:nidx
0122         plot_storage(md, idx(i), opt, 'saveit', true);
0123     end
0124 end
0125 
0126 %% extract and expand storage parameters
0127 MinStorageLevel = md.Storage.MinStorageLevel;
0128 MaxStorageLevel = md.Storage.MaxStorageLevel;
0129 if size(MinStorageLevel, 1) == 1 && ns > 1  %% expand rows
0130   MinStorageLevel = ones(ns, 1) * MinStorageLevel;
0131 end
0132 if size(MinStorageLevel, 2) == 1 && nt > 1  %% expand cols
0133   MinStorageLevel = MinStorageLevel * ones(1, nt);
0134 end
0135 if size(MaxStorageLevel, 1) == 1 && ns > 1  %% expand rows
0136   MaxStorageLevel = ones(ns, 1) * MaxStorageLevel;
0137 end
0138 if size(MaxStorageLevel, 2) == 1 && nt > 1  %% expand cols
0139   MaxStorageLevel = MaxStorageLevel * ones(1, nt);
0140 end
0141 if isempty(md.Storage.InEff)
0142   InEff = 1;                        %% no efficiency loss by default
0143 else
0144   InEff = md.Storage.InEff;
0145 end
0146 if size(InEff, 1) == 1 && ns > 1    %% expand rows
0147   InEff = ones(ns, 1) * InEff;
0148 end
0149 if size(InEff, 2) == 1 && nt > 1    %% expand cols
0150   InEff = InEff * ones(1, nt);
0151 end
0152 if isempty(md.Storage.OutEff)
0153   OutEff = 1;                       %% no efficiency loss by default
0154 else
0155   OutEff = md.Storage.OutEff;
0156 end
0157 if size(OutEff, 1) == 1 && ns > 1   %% expand rows
0158   OutEff = ones(ns, 1) * OutEff;
0159 end
0160 if size(OutEff, 2) == 1 && nt > 1   %% expand cols
0161   OutEff = OutEff * ones(1, nt);
0162 end
0163 if isempty(md.Storage.LossFactor)
0164   LossFactor = 0;                       %% no losses by default
0165 else
0166   LossFactor = md.Storage.LossFactor;
0167 end
0168 if size(LossFactor, 1) == 1 && ns > 1   %% expand rows
0169   LossFactor = ones(ns, 1) * LossFactor;
0170 end
0171 if size(LossFactor, 2) == 1 && nt > 1   %% expand cols
0172   LossFactor = LossFactor * ones(1, nt);
0173 end
0174 if isempty(md.Storage.rho)
0175   rho = 1;                      %% use worst case by default (for backward compatibility)
0176 else
0177   rho = md.Storage.rho;
0178 end
0179 if size(rho, 1) == 1 && ns > 1  %% expand rows
0180   rho = ones(ns, 1) * rho;
0181 end
0182 if size(rho, 2) == 1 && nt > 1  %% expand cols
0183   rho = rho * ones(1, nt);
0184 end
0185 
0186 %% initialize data structures to be plotted
0187 offset = opt.separation/2;          %% offset for plotting on horiz axis
0188 p = (1:nt)';
0189 pp = 2*p;
0190 ppp(pp-1) = p - offset;
0191 ppp(pp  ) = p + offset;
0192 Sp = md.results.Sp(idx, :);                     %% s+
0193 Sm = md.results.Sm(idx, :);                     %% s-
0194 eS = md.Storage.ExpectedStorageState(idx, :);   %% expected s
0195 Smin = MinStorageLevel(idx, :);     %% physical lower bound
0196 Smax = MaxStorageLevel(idx, :);     %% physical upper bound
0197 SSp   = NaN(2*nt, 1);               %% expanded for 2 pts per period
0198 SSm   = NaN(2*nt, 1);               %%    "
0199 eeS   = NaN(2*nt, 1);               %%    "
0200 SSmin(pp)   = sum(Smin, 1);
0201 SSmax(pp)   = sum(Smax, 1);
0202 SSmin(pp-1) = SSmin(pp);
0203 SSmax(pp-1) = SSmax(pp);
0204 
0205 eSi = NaN(nj_max, nt, nidx);
0206 eSf = NaN(nj_max, nt, nidx);
0207 Si_min = NaN(nj_max, nt, nidx);
0208 Si_max = NaN(nj_max, nt, nidx);
0209 dSi = zeros(nj_max, nt, nidx);  %% deltaS, total change in stored energy from injections
0210 Pg = zeros(nj_max, nt, nidx);   %% total dispatch of specified storage units
0211 jmin = zeros(1, nt);
0212 jmax = zeros(1, nt);
0213 LossCoeff = md.Delta_T * LossFactor/2;
0214 beta1 = (1-LossCoeff) ./ (1+LossCoeff);
0215 beta2 = 1 ./ (1+LossCoeff);
0216 for t = 1:nt
0217     if t == 1
0218         if md.Storage.ForceCyclicStorage
0219             SSm(t) = sum(md.Storage.InitialStorage(idx));
0220             SSp(t) = sum(md.Storage.InitialStorage(idx));
0221         else
0222             SSm(t) = sum(md.Storage.InitialStorageLowerBound(idx));
0223             SSp(t) = sum(md.Storage.InitialStorageUpperBound(idx));
0224         end
0225         eeS(t) = sum(md.Storage.InitialStorage(idx));
0226     else
0227         SSm(2*t-1) = sum(Sm(:, t-1), 1);
0228         SSp(2*t-1) = sum(Sp(:, t-1), 1);
0229         eeS(2*t-1) = sum(eS(:, t-1), 1);
0230     end
0231     %% end of period
0232     SSm(2*t)   = sum(Sm(:, t), 1);
0233     SSp(2*t)   = sum(Sp(:, t), 1);
0234     eeS(2*t)   = sum(eS(:, t), 1);
0235 
0236     vv = md.om.get_idx();
0237     for j = 1:md.idx.nj(t)
0238         dS = -md.Delta_T * ...
0239             (InEff(:,t)  .* md.QP.x(vv.i1.Psc(t,j,1)-1+(1:ns)) + ...
0240              1./OutEff(:,t) .* md.QP.x(vv.i1.Psd(t,j,1)-1+(1:ns)) );
0241         dSi(j,t,:) = dS(idx) * baseMVA;
0242         Pg(j,t,:)  = md.flow(t,j,1).mpc.gen(g, PG);
0243         Lij = md.tstep(t).Li( j:md.idx.nj(t):(ns-1)*md.idx.nj(t)+j, :);
0244         Lfj = md.tstep(t).Lf( j:md.idx.nj(t):(ns-1)*md.idx.nj(t)+j, :);
0245         Mj  = md.tstep(t).Mg( j:md.idx.nj(t):(ns-1)*md.idx.nj(t)+j, :) + ...
0246               md.tstep(t).Mh( j:md.idx.nj(t):(ns-1)*md.idx.nj(t)+j, :);
0247         Nj  = md.tstep(t).Ng( j:md.idx.nj(t):(ns-1)*md.idx.nj(t)+j, :) + ...
0248               md.tstep(t).Nh( j:md.idx.nj(t):(ns-1)*md.idx.nj(t)+j, :);
0249         temp_eSi = Lij * md.Storage.InitialStorage/baseMVA + Mj * md.QP.x;
0250         temp_eSf = Lfj * md.Storage.InitialStorage/baseMVA + Nj * md.QP.x;
0251         eSi(j,t,:) = baseMVA * temp_eSi(idx);
0252         eSf(j,t,:) = baseMVA * temp_eSf(idx);
0253         if t == 1
0254             if md.Storage.ForceCyclicStorage
0255                 Si_min(j,t,:) = md.Storage.InitialStorage(idx);
0256                 Si_max(j,t,:) = md.Storage.InitialStorage(idx);
0257             else
0258                 Si_min(j,t,:) = md.Storage.InitialStorageLowerBound(idx);
0259                 Si_max(j,t,:) = md.Storage.InitialStorageUpperBound(idx);
0260             end
0261         else
0262             tmp_eSi = reshape(eSi(j,t,:), nidx, 1);
0263             Si_min(j,t,:) = rho(idx,t) .* Sm(:,t-1) + (1-rho(idx,t)) .* tmp_eSi;
0264             Si_max(j,t,:) = rho(idx,t) .* Sp(:,t-1) + (1-rho(idx,t)) .* tmp_eSi;
0265         end
0266     end
0267 %     if nidx == -1
0268 %         rows = opt.sort_tol * round(opt.sort_tol^-1 * [ beta1(idx,t)*Si_min(:,t,1) + beta2(idx,t)*dSi(:,t,1) Si_min(:,t,1) ]);
0269 %         [junk, iii] = sortrows(rows, [1 2]);
0270 %         jmin(t) = iii(1);
0271 %         rows = opt.sort_tol * round(opt.sort_tol^-1 * [ beta1(idx,t)*Si_max(:,t,1) + beta2(idx,t)*dSi(:,t,1) Si_max(:,t,1) ]);
0272 %         [junk, iii] = sortrows(rows, [-1 -2]);
0273 %         jmax(t) = iii(1);
0274 %     else
0275         tmpS = reshape(Si_min(:,t,:), nj_max, nidx);
0276         rows = opt.sort_tol * round(opt.sort_tol^-1 * [tmpS * beta1(idx,t) + reshape(dSi(:,t,:), nj_max, nidx) * beta2(idx,t) tmpS]);
0277         [junk, i] = sortrows(rows, [1 2]);
0278         jmin(t) = i(1);
0279         tmpS = reshape(Si_max(:,t,:), nj_max, nidx);
0280         rows = opt.sort_tol * round(opt.sort_tol^-1 * [tmpS * beta1(idx,t) + reshape(dSi(:,t,:), nj_max, nidx) * beta2(idx,t) tmpS]);
0281         [junk, i] = sortrows(rows, [-1 -2]);
0282         jmax(t) = i(1);
0283 %     end
0284 end
0285 SSi_min = zeros(nt,1);
0286 SSi_max = zeros(nt,1);
0287 for t = 1:nt
0288     SSi_min(t) = sum(Si_min(jmin(t), t, :), 3);
0289     SSi_max(t) = sum(Si_max(jmax(t), t, :), 3);
0290 end
0291 
0292 %% do plots
0293 %% figure out axis limits
0294 plot(ppp, SSmax+0.001, 'LineStyle', 'none')
0295 hold on
0296 plot(ppp, SSmin-0.001, 'LineStyle', 'none')
0297 
0298 %% draw grid
0299 v = axis;
0300 if opt.show_grid
0301     m = v(3) - 100*(v(4) - v(3));
0302     M = v(4) + 100*(v(4) - v(3));
0303     for k = 0:nt
0304         line([k; k], [m; M], 'LineWidth', 0.25, 'Color', 0.9*[1 1 1]);
0305     end
0306 end
0307 axis(v);
0308 
0309 plot(ppp, SSmin, ':k', 'LineWidth', 1);
0310 plot(ppp, SSmax, ':k', 'LineWidth', 1);
0311 for t = 1:nt
0312     for j = 1:md.idx.nj(t)
0313         if opt.show_adjusted_dispatches
0314             Si = sum(Si_min(j,t,:), 3);
0315             Sf = reshape(Si_min(j,t,:), 1, nidx) * beta1(idx, t) + ...
0316                     reshape(dSi(j,t,:), 1, nidx) * beta2(idx,t);
0317             line(t+offset*[-1;1], [Si; Sf], 'Color', 0.8*[1 1 1]);
0318             Si = sum(Si_max(j,t,:), 3);
0319             Sf = reshape(Si_max(j,t,:), 1, nidx) * beta1(idx, t) + ...
0320                     reshape(dSi(j,t,:), 1, nidx) * beta2(idx,t);
0321             line(t+offset*[-1;1], [Si; Sf], 'Color', 0.8*[1 1 1]);
0322         end
0323 
0324         if opt.show_dispatches
0325             Si = sum(Si_min(j,t,:), 3);
0326             Sf = Si - sum(Pg(j,t,:), 3);
0327             line(t+offset*[-1;1], [Si; Sf], 'Color', 0.8*[1 0.7 1]);
0328             Si = sum(Si_max(j,t,:), 3);
0329             Sf = Si - sum(Pg(j,t,:), 3);
0330             line(t+offset*[-1;1], [Si; Sf], 'Color', 0.8*[1 0.7 1]);
0331         end
0332     end
0333 end
0334 plot(ppp, SSp, '--', 'LineWidth', 1, 'Color', 0.8*[0 1 0]);
0335 plot(ppp, SSm, '--r', 'LineWidth', 1, 'Color', 0.9*[1,0,0]);
0336 plot(p+offset, SSp(pp), 'v', 'LineWidth', 1, 'Color', 0.8*[0 1 0], 'MarkerSize', 6*opt.size_factor);
0337 plot(p+offset, SSm(pp), '^r', 'LineWidth', 1, 'Color', 0.9*[1,0,0], 'MarkerSize', 6*opt.size_factor);
0338 plot(ppp, eeS, 'Color', [0 0 1], 'LineWidth', 2);
0339 if opt.show_expected_initial
0340     if any(any(rho(idx,:) > 0))
0341         plot(p-offset, sum(Si_min, 3), '+', 'LineWidth', 1, 'Color', 0.9*[1,0,0], 'MarkerSize', 4*opt.size_factor);
0342         plot(p-offset, sum(Si_max, 3), '+', 'LineWidth', 1, 'Color', 0.8*[0,1,0], 'MarkerSize', 4*opt.size_factor);
0343     end
0344     plot(p-offset, SSi_min, '.', 'LineWidth', 1, 'MarkerSize', 13*opt.size_factor, 'Color', 0.9*[1 0 0]);
0345     plot(p-offset, SSi_max, '.', 'LineWidth', 1, 'MarkerSize', 13*opt.size_factor, 'Color', 0.8*[0 1 0]);
0346     plot(p-offset, sum(eSi, 3), 'o', 'LineWidth', 0.5, 'MarkerSize', 5*opt.size_factor);
0347 end
0348 if opt.show_expected_final
0349     plot(p+offset, sum(eSf, 3), 'x', 'LineWidth', 0.5, 'MarkerSize', 5*opt.size_factor);
0350 end
0351 
0352 hold off;
0353 
0354 if nidx == 1
0355     title(sprintf('Stored Energy for Storage Unit %d (Gen %d) @ Bus %d', idx, g, b), 'FontSize', 18*opt.size_factor);
0356 else
0357     if nidx == ns && all(idx == (1:ns)')
0358         txt = 'All';
0359     else
0360         txt = 'Selected';
0361     end
0362     title(sprintf('Stored Energy for %s Storage Units', txt), 'FontSize', 18*opt.size_factor);
0363 end
0364 ylabel('Energy, MWh', 'FontSize', 16*opt.size_factor);
0365 xlabel(my_xlabel, 'FontSize', 16*opt.size_factor);
0366 %legend('contract', wlabels{1}, wlabels{2}, wlabels{3}, wlabels{4}, 'Up Res', 'Dn Res', 'Lim', 'Bind', 'Location', 'EastOutside');
0367 set(gca, 'FontSize', 12*opt.size_factor);
0368 h = gcf;
0369 set(h, 'PaperOrientation', 'landscape');
0370 %set(h, 'PaperPosition', [0.25 0.25 10.5 8]);
0371 set(h, 'PaperPosition', [0 0 11 8.5]);
0372 if opt.saveit || opt.saveall
0373     if nidx == 1
0374         txt = sprintf('%d', idx);
0375     elseif nidx == ns && all(idx == (1:ns)')
0376         txt = 'all';
0377     else
0378         txt = 'selected';
0379     end
0380     if isempty(strfind(opt.savename, '%s'))
0381         pdf_name = opt.savename;
0382     else
0383         pdf_name = sprintf(opt.savename, txt);
0384     end
0385     print('-dpdf', fullfile(opt.savepath, pdf_name));
0386 end
0387 if nargout
0388     hh = h;
0389 end