MATPOWER · rdzman · Jun 4, 2024 · May 22, 2024 · May 22, 2024 · May 29, 2024
diff --git a/CHANGES.md b/CHANGES.md
@@ -2,6 +2,16 @@ Change history for MOST
 =======================
 
 
+since 1.3
+---------
+
+#### 5/29/24
+ - Fix [issue #45][12] where `most()` does not properly handle cases with
+ contingencies defined only in some periods/scenarios.
+ *Thanks to Stefano Nicolin.*
+ - Update `most_summary()` to skip display of non-existent contingencies.
+
+
 Version 1.3 - *May 10, 2024*
 ----------------------------
 
@@ -338,3 +348,4 @@ Version 1.0 - *Jun 1, 2016*
 [9]: https://github.com/MATPOWER/most/issues/37
 [10]: https://github.com/MATPOWER/most/issues/39
 [11]: https://github.com/MATPOWER/mp-opt-model
+[12]: https://github.com/MATPOWER/most/issues/45
diff --git a/docs/src/MOST-manual/MOST-manual.tex b/docs/src/MOST-manual/MOST-manual.tex
@@ -159,7 +159,7 @@
 \newcommand{\mostname}[0]{{{\bf M}{\sc atpower} \textbf{O}ptimal \textbf{S}cheduling \textbf{T}ool}}
 \newcommand{\mosturl}[0]{https://github.com/MATPOWER/most}
 \newcommand{\mostlink}[0]{\href{\mosturl}{\most{}}}
-\newcommand{\mostver}[0]{1.3}
+\newcommand{\mostver}[0]{1.3.1-dev}
 \newcommand{\md}[0]{{\most{} Data struct}}
 \newcommand{\powerweb}[0]{{\sc PowerWeb}}
 \newcommand{\pserc}[0]{{\sc PSerc}}
@@ -240,7 +240,7 @@
 %%% BEFORE PUBLISHING A NEW VERSION:
 %%% Update the publication year for \bibitem{matpower} and
 %%% \bibitem{matpower_manual} to the year of the latest release
-\date{May 10, 2024} % comment this line to display the current date
+%\date{May 10, 2024} % comment this line to display the current date
 %\date{December 14, 2011\thanks{Second revision. First revision was December 13, 2011}} % comment this line to display the current date
 
 %%% BEGIN DOCUMENT
@@ -1532,7 +1532,7 @@ \subsubsection{{\tt sd} -- Storage Data ({\tt StorageData})}
 \subsubsection{{\tt contab} -- Contingency Table}
 \label{sec:contab}
 
-The optional \code{contab} argument is a contingency table with a master set of contingencies used for security throughout entire horizon. It is a matrix in the form of a \emph{change table} recognized by \code{apply\_changes}, described in Section~\ref{MUM-sec:apply_changes} in the \mum{}. The probabilities defined in this contingency table correspond to the conditional probabilities $\psi_0^{tjk}$ of contingency~$k$ occuring conditioned on being in base scenario~$j$. While the \md{} (\code{md}) itself allows for contingencies to be defined independently for all scenarios and time periods, \code{loadmd} applies a single set of contingencies and conditional probabilities (single \code{contab}) to all.
+The optional \code{contab} argument is a contingency table with a master set of contingencies used for security throughout the entire horizon. It is a matrix in the form of a \emph{change table} recognized by \code{apply\_changes}, described in Section~\ref{MUM-sec:apply_changes} in the \mum{}. The probabilities defined in this contingency table correspond to the conditional probabilities $\psi_0^{tjk}$ of contingency~$k$ occuring conditioned on being in base scenario~$j$. While the \md{} (\code{md}) itself allows for contingencies to be defined independently for all scenarios and time periods, \code{loadmd} applies a single set of contingencies and conditional probabilities (single \code{contab}) to all.
 
 \clearpage
 
@@ -3459,10 +3459,35 @@ \subsubsection*{Incompatible Changes}
 \end{itemize}
 
 
-% \subsection{Version 1.2 -- released ??? ?, 2022}
-% \label{app:v12}
+\subsection{Version 1.3.1-dev -- released ??? ?, 202?}
+\label{app:v12}
+
+The \href{https://matpower.org/docs/MOST-manual-1.3.1.pdf}{\most{} 1.3.1 User's Manual} is available online.\footnote{\url{https://matpower.org/docs/MOST-manual-1.3.1.pdf}}
+
+\subsubsection*{Changes}
+\begin{itemize}
+\item \code{most\_summary()} now skips display of non-existent contingencies.
+\item 
+
+\end{itemize}
+
+\subsubsection*{Bugs Fixed}
+\begin{itemize}
+\item Fix issue \#45 where \code{most()} does not properly handle cases with contingencies defined only in some periods/scenarios.
+\emph{Thanks to Stefano Nicolin.}
+\item 
+\end{itemize}
+
+\subsubsection*{Incompatible Changes}
+\begin{itemize}
+\item 
+\end{itemize}
+
+
+% \subsection{Version 1.4 -- released ??? ?, 202?}
+% \label{app:v14}
 % 
-% The \href{https://matpower.org/docs/MOST-manual-1.2.pdf}{\most{} 1.2 User's Manual} is available online.\footnote{\url{https://matpower.org/docs/MOST-manual-1.2.pdf}}
+% The \href{https://matpower.org/docs/MOST-manual-1.4.pdf}{\most{} 1.4 User's Manual} is available online.\footnote{\url{https://matpower.org/docs/MOST-manual-1.4.pdf}}
 % 
 % \subsubsection*{Changes}
 % \begin{itemize}
@@ -3495,7 +3520,7 @@ \subsubsection*{Incompatible Changes}
 \doi{10.1109/TPWRS.2010.2051168}
 
 \bibitem{matpower}
-R.~D. Zimmerman, C.~E. Murillo-S{\'a}nchez (2022). \matpower{}\\~
+R.~D. Zimmerman, C.~E. Murillo-S{\'a}nchez (2024). \matpower{}\\~
 [Software]. Available: \url{https://matpower.org}\\
 \doi{10.5281/zenodo.3236535}
 

diff --git a/lib/most.m b/lib/most.m
@@ -144,9 +144,19 @@
  mo.IncludeFixedReserves = 0;
  end
 end
+if mo.SecurityConstrained % check that at least some contingency is specified
+ have_contingency = 0;
+ for t = 1:nt
+ for j = 1:mdi.idx.nj(t)
+ if isfield(mdi, 'cont') && isfield(mdi.cont(t,j), 'contab') && ...
+ ~isempty(mdi.cont(t,j).contab)
+ have_contingency = 1; % found a contingency
+ end
+ end
+ end
+end
 if mo.SecurityConstrained == -1
- if isfield(mdi, 'cont') && isfield(mdi.cont(1,1), 'contab') && ...
- ~isempty(mdi.cont(1,1).contab)
+ if have_contingency
  mo.SecurityConstrained = 1;
  else
  mo.SecurityConstrained = 0;
@@ -165,9 +175,8 @@
  isfield(mdi.FixedReserves(1,1,1), 'req'))
  error('most: MDI.FixedReserves(t,j,k) must be specified when MPOPT.most.fixed_res = 1');
 end
-if mo.SecurityConstrained && ~(isfield(mdi, 'cont') && ...
- isfield(mdi.cont(1,1), 'contab') && ~isempty(mdi.cont(1,1).contab))
- error('most: MDI.cont(t,j).contab cannot be empty when MPOPT.most.security_constraints = 1');
+if mo.SecurityConstrained && ~have_contingency
+ error('most: MDI.cont(t,j).contab cannot be empty for all t, j when MPOPT.most.security_constraints = 1');
 end
 if mo.IncludeFixedReserves && mo.SecurityConstrained
  warning('most: Using MPOPT.most.fixed_res = 1 and MPOPT.most.security_constraints = 1 together is not recommended.');
@@ -446,35 +455,40 @@
  mdi.StepProb(t) = sum(scenario_probs); % probability of making it to the t-th step
  if mdi.SecurityConstrained
  for j = 1:mdi.idx.nj(t)
- [tmp, ii] = sort(mdi.cont(t,j).contab(:, CT_LABEL)); %sort in ascending contingency label
- contab = mdi.cont(t,j).contab(ii, :);
- rowdecomlist = ones(size(contab,1), 1);
- for l = 1:size(contab, 1)
- if contab(l, CT_TABLE) == CT_TGEN && contab(l, CT_COL) == GEN_STATUS ...
- && contab(l, CT_CHGTYPE) == CT_REP && contab(l, CT_NEWVAL) == 0 ... % gen turned off
- && mdi.flow(t,j,1).mpc.gen(contab(l, CT_ROW), GEN_STATUS) <= 0 % but it was off on input
- rowdecomlist(l) = 0;
- elseif contab(l, CT_TABLE) == CT_TBRCH && contab(l, CT_COL) == BR_STATUS ...
- && contab(l, CT_CHGTYPE) == CT_REP && contab(l, CT_NEWVAL) == 0 ... % branch taken out
- && mdi.flow(t,j,1).mpc.branch(contab(l, CT_ROW), BR_STATUS) <= 0 % but it was off on input
- rowdecomlist(l) = 0;
+ if isempty(mdi.cont(t,j).contab)
+ mdi.idx.nc(t, j) = 0;
+ mdi.CostWeights(1, j, t) = 1;
+ else
+ [tmp, ii] = sort(mdi.cont(t,j).contab(:, CT_LABEL)); %sort in ascending contingency label
+ contab = mdi.cont(t,j).contab(ii, :);
+ rowdecomlist = ones(size(contab,1), 1);
+ for l = 1:size(contab, 1)
+ if contab(l, CT_TABLE) == CT_TGEN && contab(l, CT_COL) == GEN_STATUS ...
+ && contab(l, CT_CHGTYPE) == CT_REP && contab(l, CT_NEWVAL) == 0 ... % gen turned off
+ && mdi.flow(t,j,1).mpc.gen(contab(l, CT_ROW), GEN_STATUS) <= 0 % but it was off on input
+ rowdecomlist(l) = 0;
+ elseif contab(l, CT_TABLE) == CT_TBRCH && contab(l, CT_COL) == BR_STATUS ...
+ && contab(l, CT_CHGTYPE) == CT_REP && contab(l, CT_NEWVAL) == 0 ... % branch taken out
+ && mdi.flow(t,j,1).mpc.branch(contab(l, CT_ROW), BR_STATUS) <= 0 % but it was off on input
+ rowdecomlist(l) = 0;
+ end
  end
+ contab = contab(rowdecomlist ~= 0, :);
+ mdi.cont(t, j).contab = contab;
+ clist = unique(contab(:, CT_LABEL));
+ mdi.idx.nc(t, j) = length(clist);
+ k = 2;
+ for label = clist'
+ mdi.flow(t, j, k).mpc = apply_changes(label, mdi.flow(t, j, 1).mpc, contab);
+ ii = find( label == contab(:, CT_LABEL) );
+ mdi.CostWeights(k, j, t) = contab(ii(1), CT_PROB);
+ mdi.idx.nb(t, j, k) = size(mdi.flow(t, j, k).mpc.bus, 1);
+ mdi.idx.ny(t, j, k) = length(find(mdi.flow(t, j, 1).mpc.gencost(:, MODEL) == PW_LINEAR));
+ k = k + 1;
+ end
+ mdi.CostWeights(1, j, t) = 1 - sum(mdi.CostWeights(2:mdi.idx.nc(t,j)+1, j, t));
+ mdi.CostWeights(1:mdi.idx.nc(t,j)+1, j, t) = scenario_probs(j) * mdi.CostWeights(1:mdi.idx.nc(t,j)+1, j, t);
  end
- contab = contab(rowdecomlist ~= 0, :);
- mdi.cont(t, j).contab = contab;
- clist = unique(contab(:, CT_LABEL));
- mdi.idx.nc(t, j) = length(clist);
- k = 2;
- for label = clist'
- mdi.flow(t, j, k).mpc = apply_changes(label, mdi.flow(t, j, 1).mpc, contab);
- ii = find( label == contab(:, CT_LABEL) );
- mdi.CostWeights(k, j, t) = contab(ii(1), CT_PROB);
- mdi.idx.nb(t, j, k) = size(mdi.flow(t, j, k).mpc.bus, 1);
- mdi.idx.ny(t, j, k) = length(find(mdi.flow(t, j, 1).mpc.gencost(:, MODEL) == PW_LINEAR));
- k = k + 1;
- end
- mdi.CostWeights(1, j, t) = 1 - sum(mdi.CostWeights(2:mdi.idx.nc(t,j)+1, j, t));
- mdi.CostWeights(1:mdi.idx.nc(t,j)+1, j, t) = scenario_probs(j) * mdi.CostWeights(1:mdi.idx.nc(t,j)+1, j, t);
  end
  else
  for j = 1:mdi.idx.nj(t)

diff --git a/lib/most_summary.m b/lib/most_summary.m
@@ -68,28 +68,30 @@
 nl = size(mpc.branch, 1);
 ng = size(mpc.gen, 1);
 nt = mdo.idx.nt;
-nj_max = max(mdo.idx.nj);
-nc_max = max(max(mdo.idx.nc));
+nj = mdo.idx.nj;
+nc = mdo.idx.nc;
+nj_max = max(nj);
+nc_max = max(max(nc));
 ns = mdo.idx.ns;
 
 %% summarize results
-psi = zeros(nt, nj_max, nc_max+1);
-Pg = zeros(ng, nt, nj_max, nc_max+1);
-Pd = zeros(nb, nt, nj_max, nc_max+1);
+psi = NaN(nt, nj_max, nc_max+1);
+Pg = NaN(ng, nt, nj_max, nc_max+1);
+Pd = NaN(nb, nt, nj_max, nc_max+1);
 if mdo.idx.ntramp
  Rup = mdo.results.Rrp;
  Rdn = mdo.results.Rrm;
 else
  Rup = [];
  Rdn = [];
 end
-u = zeros(ng, nt);
-lamP = zeros(nb, nt, nj_max, nc_max+1);
-muF = zeros(nl, nt, nj_max, nc_max+1);
-Pf = zeros(nl, nt, nj_max, nc_max+1);
+u = NaN(ng, nt);
+lamP = NaN(nb, nt, nj_max, nc_max+1);
+muF = NaN(nl, nt, nj_max, nc_max+1);
+Pf = NaN(nl, nt, nj_max, nc_max+1);
 for t = 1:nt
- for j = 1:mdo.idx.nj(t)
- for k = 1:mdo.idx.nc(t,j)+1
+ for j = 1:nj(t)
+ for k = 1:nc(t,j)+1
  rr = mdo.flow(t,j,k).mpc;
  psi(t, j, k) = mdo.CostWeightsAdj(k, j, t);
  u(:, t) = rr.gen(:, GEN_STATUS);
@@ -159,19 +161,19 @@
  end
  fprintf('\n');
 
- print_most_summary_section('PG', 'Gen', nt, nj_max, nc_max, Pg);
+ print_most_summary_section('PG', 'Gen', nt, nj_max, nc, Pg);
  if mdo.idx.ntramp
  print_most_summary_section('RAMP UP', 'Gen', nt, 1, 0, Rup);
  print_most_summary_section('RAMP DOWN', 'Gen', nt, 1, 0, Rdn);
  end
- print_most_summary_section('FIXED LOAD', 'Bus', nt, nj_max, nc_max, Pd);
+ print_most_summary_section('FIXED LOAD', 'Bus', nt, nj_max, nc, Pd);
  if ns
  print_most_summary_section('ESS E[SoC]', 'ESS', nt, 1, 0, SoC);
  end
  if mdo.DCMODEL
- print_most_summary_section('LAM_P', 'Bus', nt, nj_max, nc_max, lamP);
- print_most_summary_section('PF', 'Brch', nt, nj_max, nc_max, Pf);
- print_most_summary_section('MU_F', 'Brch', nt, nj_max, nc_max, muF);
+ print_most_summary_section('LAM_P', 'Bus', nt, nj_max, nc, lamP);
+ print_most_summary_section('PF', 'Brch', nt, nj_max, nc, Pf);
+ print_most_summary_section('MU_F', 'Brch', nt, nj_max, nc, muF);
  end
 end
 
@@ -180,7 +182,7 @@
 end
 
 %%---------------------------------------------------------
-function print_most_summary_section(label, section_type, nt, nj_max, nc_max, data, tol)
+function print_most_summary_section(label, section_type, nt, nj_max, nc, data, tol)
 if nargin < 7
  tol = 1e-4;
 end
@@ -189,6 +191,7 @@ function print_most_summary_section(label, section_type, nt, nj_max, nc_max, dat
 fprintf('\n==========%-12s==========\n', sprintf('%s%s', bl, label));
 if any(data(:))
  for j = 1:nj_max
+ nc_max = max(nc(:,j));
  for k = 1:nc_max+1
  if nj_max > 1 || nc_max > 0
  fprintf('\nSCENARIO %d', j);

diff --git a/lib/mostver.m b/lib/mostver.m
@@ -24,9 +24,9 @@
 % See https://github.com/MATPOWER/most for more info.
 
 v = struct( 'Name', 'MOST', ...
- 'Version', '1.3', ...
+ 'Version', '1.3.1-dev', ...
  'Release', '', ...
- 'Date', '10-May-2024' );
+ 'Date', '29-May-2024' );
 if nargout > 0
  if nargin > 0
  rv = v;