git subrepo pull mp-opt-model

subrepo: subdir: "mp-opt-model" merged: "03a6be78" upstream: origin: "[email protected]:MATPOWER/mp-opt-model" branch: "master" commit: "03a6be78" git-subrepo: version: "0.4.6" origin: "https://github.com/ingydotnet/git-subrepo" commit: "110b9eb"
MATPOWER · Feb 1, 2024 · fdd7dfe · fdd7dfe
1 parent 4f22b4f
commit fdd7dfe
Show file tree

Hide file tree

Showing 35 changed files with 815 additions and 628 deletions.
diff --git a/mp-opt-model/.gitrepo b/mp-opt-model/.gitrepo
@@ -6,7 +6,7 @@
 [subrepo]
 	remote = [email protected]:MATPOWER/mp-opt-model
 	branch = master
-	commit = 478c42c9b74f8dc6dabeb5a073e4a1aa0196d7e8
-	parent = 07237f516bf7550a1df7b2c9ca3a1f63accb2078
+	commit = 03a6be782f3c0ba691e56a8d2da255a6d1974708
+	parent = 4f22b4f04b042c20cc78672d578d0cd28f944a5e
 	method = merge
-	cmdver = 0.4.5
+	cmdver = 0.4.6
diff --git a/mp-opt-model/CHANGES.md b/mp-opt-model/CHANGES.md
@@ -5,6 +5,48 @@ Change history for MP-Opt-Model
 Since version 4.1
 -----------------
 
+#### 1/31/24
+  - Add `convert_lin_constraint()` and `convert_lin_constraint_multipliers()`
+    functions to eliminate code duplication for common task of converting
+    linear constraints and their multipliers between a single set of
+    doubly-bounded inequality constraints and separate sets of equality
+    and upper-bounded inequality constraints.
+  - Change solver for CPLEX price computation stage in `miqps_cplex()` from
+    primal simplex to dual simplex (probably no impact, except it was a
+    simple way to get a newly failing test in another project to pass again).
+
+#### 12/8/23
+  - Always skip price computation stage in `miqps_<solver>()` functions for
+    pure (as opposed to mixed) integer problems.
+
+#### 11/29/23
+  - Add support to `miqps_master()` for calling `miqps_<my_solver>()` by
+    setting `opt.alg` to `'<MY_SOLVER>'` to allow for handling custom
+    MILP/MIQP solvers.
+
+#### 11/8/23
+  - Add support to `nlps_master()` for calling `nlps_<my_solver>()` by setting
+    `opt.alg` to `'<MY_SOLVER>'` to allow for handling custom NLP solvers.
+
+#### 11/6/23
+  - Add to `opt_model/add_lin_constraint()` the option to provide/store
+    the transpose of the `A` matrix instead of the original. This can
+    potentially save significant memory for sparse matrices with many more
+    columns than rows. E.g. storage constraints in [MOST][11] for 8760 hour
+    planning horizon.
+
+#### 10/13/23
+  - Update `opt_model/display_soln()` to avoid displaying an infinite
+    average for quadratic costs when corresponding quantity is zero.
+
+#### 9/13/23
+  - Clear cached parameters after updating quadratic costs via
+    `opt_model/set_params()`.
+
+#### 9/12/23
+  - Tweak some MI/QPS solver tests to make them more robust for `'DEFAULT'`
+    solver under different environments.
+
 #### 3/27/23
   - Update for compatibility with MATLAB R2023a (Optimization Toolbox 9.5)
     which removed `x0` as a valid input to `linprog`.
@@ -32,7 +74,7 @@ Version 4.1 - *Dec 13, 2022*
   - Add `runtime` field to `output` argument of `qps_glpk()` and
     `qps_mosek()`.
   - Add support to `qps_master()` for calling `qps_<my_solver>()` by setting
-    `opt.alg` to `'<MY_SOLVER>'` to allow for custom solvers.
+    `opt.alg` to `'<MY_SOLVER>'` to allow for handling custom LP/QP solvers.
 
 #### 7/5/22
   - Update for compatibility with Artelys Knitro 13.1 and later.
@@ -425,3 +467,4 @@ Version 0.7.0 - *Jun 20, 2019*
 [8]: https://github.com/MATPOWER/mptest
 [9]: https://github.com/MATPOWER/mips
 [10]: https://savannah.gnu.org/bugs/?52614
+[11]: https://github.com/MATPOWER/most
diff --git a/mp-opt-model/docs/src/MP-Opt-Model-manual/MP-Opt-Model-manual.tex b/mp-opt-model/docs/src/MP-Opt-Model-manual/MP-Opt-Model-manual.tex
@@ -151,10 +151,11 @@
 \newcommand{\mpomlink}[0]{\href{\mpomurl}{\mpom{}}}
 \newcommand{\mpomname}[0]{\mpom{}}
 % \newcommand{\mpomname}[0]{{{\bf M}{\sc at}{\bf P}{\sc ower} \textbf{Opt}imization \textbf{Model}}}
-\newcommand{\mpomver}[0]{4.1}
+\newcommand{\mpomver}[0]{4.1+}
 \newcommand{\most}[0]{{MOST}}
 \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.2}
 \newcommand{\syngrid}[0]{{SynGrid}}
 \newcommand{\syngridver}[0]{1.0.1}
@@ -242,7 +243,7 @@
 %\title{\hl{--- DRAFT  ---}\\\hl{\em do not distribute}\\~\\{\huge \bfseries \mpomname{} User's Manual } \\ ~ \\ \LARGE Version \mpomver{}\\
 \title{{\huge \bfseries \mpomname{} User's Manual } \\ ~ \\ \LARGE Version \mpomver{}}
 \author{Ray~D.~Zimmerman}
-\date{December 13, 2022} % comment this line to display the current date
+%\date{December 13, 2022} % comment this line to display the current date
 %\date{May 8, 2020\thanks{Second revision. First revision was April 29, 2020}} % comment this line to display the current date
 
 %%% BEGIN DOCUMENT
@@ -253,7 +254,7 @@
 \vfill
 \begin{center}
 {\scriptsize
-\copyright~2020--2022~\PSERC{}\\
+\copyright~2020--2023~\PSERC{}\\
 All Rights Reserved}
 \end{center}
 
@@ -1026,7 +1027,7 @@ \subsection{NLP Solvers -- {\tt nlps\_master}}
 \end{threeparttable}
 \end{table}
 
-The user-defined functions for evaluating the objective function, constraints and Hessian are identical to those required by \mipslink{}j. That is, they identical to those required by \code{fmincon}, with one exception described below for the Hessian evaluation function. Specifically, \code{f\_fcn} should return \code{f} as the scalar objective function value $f(x)$, \code{df} as an $n \times 1$ vector equal to $\nabla f$ and, unless \code{gh\_fcn} is provided and the Hessian is computed by \code{hess\_fcn}, \code{d2f} as an $n \times n$ matrix equal to the Hessian $\der{^2f}{x^2}$. Similarly, the constraint evaluation function \code{gh\_fcn} must return the $m \times 1$ vector of nonlinear equality constraint violations $g(x)$, the $p \times 1$ vector of nonlinear inequality constraint violations $h(x)$ along with their gradients in \code{dg} and \code{dh}. Here \code{dg} is an $n \times m$ matrix whose $j^\mathrm{th}$ column is $\nabla g_j$ and \code{dh} is $n \times p$, with $j^\mathrm{th}$ column equal to $\nabla h_j$. Finally, for cases with nonlinear constraints, \code{hess\_fcn} returns the $n \times n$ Hessian $\der{^2\mathcal{L}}{x^2}$ of the Lagrangian function
+The user-defined functions for evaluating the objective function, constraints and Hessian are identical to those required by \mipslink{}. That is, they are identical to those required by \code{fmincon}, with one exception described below for the Hessian evaluation function. Specifically, \code{f\_fcn} should return \code{f} as the scalar objective function value $f(x)$, \code{df} as an $n \times 1$ vector equal to $\nabla f$ and, unless \code{gh\_fcn} is provided and the Hessian is computed by \code{hess\_fcn}, \code{d2f} as an $n \times n$ matrix equal to the Hessian $\der{^2f}{x^2}$. Similarly, the constraint evaluation function \code{gh\_fcn} must return the $m \times 1$ vector of nonlinear equality constraint violations $g(x)$, the $p \times 1$ vector of nonlinear inequality constraint violations $h(x)$ along with their gradients in \code{dg} and \code{dh}. Here \code{dg} is an $n \times m$ matrix whose $j^\mathrm{th}$ column is $\nabla g_j$ and \code{dh} is $n \times p$, with $j^\mathrm{th}$ column equal to $\nabla h_j$. Finally, for cases with nonlinear constraints, \code{hess\_fcn} returns the $n \times n$ Hessian $\der{^2\mathcal{L}}{x^2}$ of the Lagrangian function
 \begin{equation}
 \mathcal{L}(x, \lambda, \mu, \sigma) = \sigma f(x) + \trans{\lambda} g(x) + \trans{\mu} h(x)
 \end{equation}
@@ -2382,7 +2383,7 @@ \subsubsection{General Nonlinear Costs}
 [f, df] = fcn(x)
 [f, df, d2f] = fcn(x)
 \end{Code}
-where \code{f} is a scalar with the value of the function $f(x)$, \code{df} is the $1 \times n_x$ gradient of $f$, and \code{d2f} is the $n_x \times n_x$ Hessian $H$, where $n_x$ is the number of elements in $x$.
+where \code{f} is a scalar with the value of the function $f(x)$, \code{df} is the $n_x \times 1$ gradient of $f$, and \code{d2f} is the $n_x \times n_x$ Hessian $H$, where $n_x$ is the number of elements in $x$.
 
 The first input argument \code{x} takes one of two forms. If the constraint set is added with \code{varsets} empty or missing, then \code{x} will be the full optimization vector. Otherwise it will be a cell array of vectors corresponding to the variable sets specified in \code{varsets}.
 
@@ -3169,7 +3170,7 @@ \subsubsection*{Examples:}
 \subsection{Indexed Sets}
 \label{sec:indexed_sets}
 
-A variable, constraint or cost set is typically identified simply by a \code{name}, but it is also possible to used indexed names. For example, an optimal scheduling problem with a one week horizon might include a vector variable \textbf{y} for each day, indexed from 1 to 7, and another vector variable \textbf{z} for each hour of each day, indexed from (1, 1) to (7, 24).
+A variable, constraint or cost set is typically identified simply by a \code{name}, but it is also possible to use indexed names. For example, an optimal scheduling problem with a one week horizon might include a vector variable \textbf{y} for each day, indexed from 1 to 7, and another vector variable \textbf{z} for each hour of each day, indexed from (1, 1) to (7, 24).
 
 In this case, we case use a single indexed named set for \textbf{y} and another for \textbf{z}. The dimensions are initialized via the \code{init\_indexed\_name} method before adding the variables to the model.\footnote{The same is true for indexed named sets of constraints or costs.}
 
@@ -3180,8 +3181,8 @@ \subsubsection*{\code{init\_indexed\_name}}
 
 \noindent Examples:
 \begin{Code}
-[f, df, d2f] = om.init_indexed_name('var', 'y', {7});
-[f, df, d2f] = om.init_indexed_name('var', 'z', {7, 24});
+om.init_indexed_name('var', 'y', {7});
+om.init_indexed_name('var', 'z', {7, 24});
 \end{Code}
 
 After initializing the dimensions, indexed named sets of variables, constraints or costs can be added by supplying the indices in the \code{idx\_list} argument following the \code{name} argument in the call to the corresponding \code{add\_var}, \code{add\_lin\_constraint}, \code{add\_nln\_constraint}, \code{add\_quad\_cost}, or \code{add\_nln\_cost} method. The \code{idx\_list} argument is simply a cell array containing the indices of interest.
@@ -3517,6 +3518,43 @@ \subsubsection{Methods}
 \clearpage
 \section{Utility Functions}
 
+\subsection{\tt convert\_lin\_constraint}
+\begin{Code}
+  [ieq, igt, ilt, Ae, be, Ai, bi] = convert_lin_constraint(A, l, u)
+  [ieq, igt, ilt, A, b] = convert_lin_constraint(A, l, u)
+\end{Code}
+
+This function converts linear constraints from a single set of doubly-bounded inequality constraints
+\begin{equation}
+l \le A x \le u
+\end{equation}
+to separate sets of equality and upper-bounded inequality constraints.
+\begin{align}
+A_e x &= b_e \\
+A_i x &\le b_i
+\end{align}
+
+The first three return values are index vectors which satisfy the following.
+\begin{Code}
+Ae = A(ieq, :);
+be = u(ieq, 1);
+Ai  = [ A(ilt, :); -A(igt, :) ];
+bi  = [ u(ilt, 1); -l(igt, 1) ];
+\end{Code}
+Alternatively, the returned matrices and right hand side vectors can be stacked into a single set with the equalities first, then the inequalities.
+\begin{Code}
+A = [Ae; Ai]
+b = [be; bi]
+\end{Code}
+
+\subsection{\tt convert\_lin\_constraint\_multipliers}
+\begin{Code}
+  [mu_l, mu_u] = convert_lin_constraint_multipliers(lam, mu, ieq, igt, ilt)
+\end{Code}
+
+This function converts the multipliers on linear constraints from separate sets for equality and upper-bounded inequality constraints to those for doubly-bounded inequality constraints.
+
+
 \subsection{\tt have\_fcn}
 
 This function is deprecated. Instead, please use \code{have\_feature}, now included as part of \mptestlink{} and described in the \mptest{} \href{\mptesturl/blob/master/README.md}{\code{README}} file. It is simply a drop-in replacement that has been reimplemented with an extensible, modular design, where the detection of a feature named \code{<tag>} is implemented by the function named \code{have\_feature\_<tag>}. The current \code{have\_fcn} is a simple wrapper around \code{have\_feature}.
@@ -4474,7 +4512,7 @@ \subsection{Version 4.1 -- released Dec 13, 2022}
 
 \subsubsection*{New Features}
 \begin{itemize}
-\item Add support to \code{qps\_master()} for calling \code{qps\_<my\_solver>()} by setting \code{opt.alg} to \codeq{<MY\_SOLVER>} to allow for custom solvers.
+\item Add support to \code{qps\_master()} for calling \code{qps\_<my\_solver>()} by setting \code{opt.alg} to \codeq{<MY\_SOLVER>} to allow for handling custom LP/QP solvers.
 \end{itemize}
 
 \subsubsection*{Other Changes}
@@ -4492,10 +4530,13 @@ \subsection{Version 4.2 -- released ??? ?, 202?}
 
 \subsubsection*{New Features}
 \begin{itemize}
-\item 
+\item Option for \code{opt\_model/add\_lin\_constraint()} to provide/store the transpose of the $A$ matrix instead of the original. This can potentially save significant memory for sparse matrices with many more columns than rows. E.g. storage constraints in \mostlink{} for 8760 hour planning horizon.
+\item Add support to \code{nlps\_master()} for calling \code{nlps\_<my\_solver>()} by setting \code{opt.alg} to \codeq{<MY\_SOLVER>} to allow for handling custom NLP solvers.
+\item Add support to \code{miqps\_master()} for calling \code{miqps\_<my\_solver>()} by setting \code{opt.alg} to \codeq{<MY\_SOLVER>} to allow for handling custom MILP/MIQP solvers.
 \item New functions:
     \begin{itemize}
-    \item \code{foobar()} does whizbang.
+    \item \code{convert\_lin\_constraint()} converts linear constraints from a single set of doubly-bounded inequality constraints to separate sets of equality and upper-bounded inequality constraints.
+    \item \code{convert\_lin\_constraint\_multipliers()} converts multipliers on linear constraints from separate sets for equality and upper-bounded inequality constraints to those for doubly-bounded inequality constraints.
     \end{itemize}
 \item New \code{opt\_model} methods:
     \begin{itemize}
@@ -4506,13 +4547,16 @@ \subsubsection*{New Features}
 
 \subsubsection*{Bugs Fixed}
 \begin{itemize}
+\item Clear cached parameters after updating quadratic costs via \code{opt\_model/set\_params()}.
 \item 
     \emph{Thanks to Fulano.}
 \end{itemize}
 
 \subsubsection*{Other Changes}
 \begin{itemize}
 \item Update for compatibility with MATLAB R2023a (Optimization Toolbox 9.5) which removed \code{x0} as a valid input to \code{linprog}.
+\item Always skip price computation stage in \code{miqps\_<solver>()} functions for pure (as opposed to mixed) integer problems.
+
 
 \end{itemize}
 

diff --git a/mp-opt-model/lib/@mp_idx_manager/describe_idx.m b/mp-opt-model/lib/@mp_idx_manager/describe_idx.m
@@ -23,7 +23,9 @@
 
 label = cell(size(idxs));       %% pre-allocate return cell array
 
-s = obj.set_type_idx_map(set_type, idxs);
+if ~isempty(idxs)
+    s = obj.set_type_idx_map(set_type, idxs);
+end
 for i = 1:length(idxs(:))
     idx = s(i).idx;
     if isempty(idx)

diff --git a/mp-opt-model/lib/@opt_model/add_lin_constraint.m b/mp-opt-model/lib/@opt_model/add_lin_constraint.m
@@ -1,19 +1,23 @@
-function om = add_lin_constraint(om, name, idx, A, l, u, varsets)
+function om = add_lin_constraint(om, name, idx, A, l, u, varsets, tr)
 %ADD_LIN_CONSTRAINT  Adds a set of linear constraints to the model.
 %
 %   OM.ADD_LIN_CONSTRAINT(NAME, A, L, U);
 %   OM.ADD_LIN_CONSTRAINT(NAME, A, L, U, VARSETS);
+%   OM.ADD_LIN_CONSTRAINT(NAME, A, L, U, VARSETS, TR);
 %
 %   OM.ADD_LIN_CONSTRAINT(NAME, IDX_LIST, A, L, U);
 %   OM.ADD_LIN_CONSTRAINT(NAME, IDX_LIST, A, L, U, VARSETS);
+%   OM.ADD_LIN_CONSTRAINT(NAME, IDX_LIST, A, L, U, VARSETS, TR);
 %
 %   Linear constraints are of the form L <= A * x <= U, where x is a
 %   vector made of the vars specified in VARSETS (in the order given).
 %   This allows the A matrix to be defined only in terms of the relevant
 %   variables without the need to manually create a lot of zero columns.
 %   If VARSETS is empty, x is taken to be the full vector of all
 %   optimization variables. If L or U are empty, they are assumed to be
-%   appropriately sized vectors of -Inf and Inf, respectively.
+%   appropriately sized vectors of -Inf and Inf, respectively. If TR is
+%   present and true, it means that A' is supplied/stored rather than A.
+%   In some contexts this can be used to save significant memory.
 %
 %   Indexed Named Sets
 %       A constraint set can be identified by a single NAME, as described
@@ -44,7 +48,7 @@
 %   See also OPT_MODEL, PARAMS_LIN_CONSTRAINT.
 
 %   MP-Opt-Model
-%   Copyright (c) 2008-2020, Power Systems Engineering Research Center (PSERC)
+%   Copyright (c) 2008-2023, Power Systems Engineering Research Center (PSERC)
 %   by Ray Zimmerman, PSERC Cornell
 %
 %   This file is part of MP-Opt-Model.
@@ -53,10 +57,18 @@
 
 %% initialize input arguments
 if iscell(idx)          %% indexed named set
-    if nargin < 7
-        varsets = {};
+    if nargin < 8
+        tr = 0;
+        if nargin < 7
+            varsets = {};
+        end
     end
 else                    %% simple named set
+    if nargin < 7
+        tr = 0;
+    else
+        tr = varsets;
+    end
     if nargin < 6
         varsets = {};
     else
@@ -68,7 +80,11 @@
     idx = {};
 end
 
-[N, M] = size(A);
+if tr
+    [M, N] = size(A);
+else
+    [N, M] = size(A);
+end
 if isempty(l)                   %% default l is -Inf
     l = -Inf(N, 1);
 elseif N > 1 && length(l) == 1  %% expand from scalar as needed
@@ -95,4 +111,4 @@
 end
 
 %% add the named linear constraint set
-om.add_named_set('lin', name, idx, N, A, l, u, varsets);
+om.add_named_set('lin', name, idx, N, A, l, u, varsets, tr);
diff --git a/mp-opt-model/lib/@opt_model/add_named_set.m b/mp-opt-model/lib/@opt_model/add_named_set.m
@@ -11,7 +11,7 @@
 %       OM.ADD_NAMED_SET('var', NAME, IDX_LIST, N, V0, VL, VU, VT);
 %
 %   Linear Constraint Set
-%       OM.ADD_NAMED_SET('lin', NAME, IDX_LIST, N, A, L, U, VARSETS);
+%       OM.ADD_NAMED_SET('lin', NAME, IDX_LIST, N, A, L, U, VARSETS, TR);
 %
 %   Nonlinear Inequality Constraint Set
 %       OM.ADD_NAMED_SET('nle', NAME, IDX_LIST, N, FCN, HESS, COMPUTED_BY, VARSETS);
@@ -26,7 +26,7 @@
 %            ADD_QUAD_COST and ADD_NLN_COST.
 
 %   MP-Opt-Model
-%   Copyright (c) 2008-2020, Power Systems Engineering Research Center (PSERC)
+%   Copyright (c) 2008-2023, Power Systems Engineering Research Center (PSERC)
 %   by Ray Zimmerman, PSERC Cornell
 %
 %   This file is part of MP-Opt-Model.
@@ -61,16 +61,18 @@
             om_ff.data.vt = subsasgn(om_ff.data.vt, sc, vt);    %% variable type
         end
     case 'lin'          %% linear constraint set
-        [A, l, u, varsets] = deal(varargin{:});
+        [A, l, u, varsets, tr] = deal(varargin{:});
         if isempty(idx)
             om_ff.data.A.(name)  = A;
             om_ff.data.l.(name)  = l;
             om_ff.data.u.(name)  = u;
+            om_ff.data.tr.(name)  = tr;
             om_ff.data.vs.(name) = varsets;
         else
             om_ff.data.A  = subsasgn(om_ff.data.A, sc, A);
             om_ff.data.l  = subsasgn(om_ff.data.l, sc, l);
             om_ff.data.u  = subsasgn(om_ff.data.u, sc, u);
+            om_ff.data.tr  = subsasgn(om_ff.data.tr, sc, tr);
             om_ff.data.vs = subsasgn(om_ff.data.vs, sc, varsets);
         end
         if ~isempty(om_ff.params)       %% clear cache of aggregated params

diff --git a/mp-opt-model/lib/@opt_model/add_nln_cost.m b/mp-opt-model/lib/@opt_model/add_nln_cost.m
@@ -18,7 +18,7 @@
 %       F = FCN(X)
 %       [F, DF] = FCN(X)
 %       [F, DF, D2F] = FCN(X)
-%   where F is a scalar with the value of the function, DF is the 1 x NX
+%   where F is a scalar with the value of the function, DF is the NX x 1
 %   gradient, and D2F is the NX x NX Hessian and NX is the number of
 %   elements in X.
 %

diff --git a/mp-opt-model/lib/@opt_model/display_soln.m b/mp-opt-model/lib/@opt_model/display_soln.m
@@ -141,10 +141,15 @@
                     if len == 1
                         c_constant = kk;
                         c_linear = cc' * xx;
-                        c_average = c_total / sum(xx);
+                        if abs(sum(xx)) > 1e-9
+                            c_average = c_total / sum(xx);
+                        else
+                            c_average = NaN;
+                        end
                     else
                         c_linear = cc .* xx;
                         c_average = c_total ./ xx;
+                        c_average(isinf(c_average)) = NaN;
                         if isscalar(kk)
                             c_constant = ones(len, 1)*kk/len;
                         else