diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index bf13637b8..3a49d0969 100644 --- a/.github/workflows/test_pr_and_main.yml +++ b/.github/workflows/test_pr_and_main.yml @@ -1262,7 +1262,9 @@ jobs: mpisppy/tests/test_prox_approx.py \ mpisppy/tests/test_sep_rho.py \ mpisppy/tests/test_reduced_costs_fixer.py \ - mpisppy/tests/test_slammer.py + mpisppy/tests/test_slammer.py \ + mpisppy/tests/test_dual_certificate.py \ + mpisppy/tests/test_ipopt_outer_bound.py - name: Upload coverage data if: always() @@ -1273,6 +1275,128 @@ jobs: if-no-files-found: ignore include-hidden-files: true + ipopt-tests: + name: Ipopt tests (idaes-ext build, with HSL) + runs-on: ubuntu-latest + needs: [ruff] + steps: + - uses: actions/checkout@v3 + - uses: conda-incubator/setup-miniconda@v3 + with: + activate-environment: test_env + python-version: 3.11 + - name: Install dependencies + run: | + conda install mpi4py pandas setuptools + + - name: Install Ipopt's system libraries + # The idaes-ext binary carries no RPATH and resolves libgfortran, + # liblapack and libblas from the system. A bare runner has none of + # them, and the failure surfaces as an unhelpful load error at first + # solve rather than at install. Pyomo's workflow installs the same + # three for the same reason. + run: | + sudo apt-get update + sudo apt-get install -y libopenblas-dev gfortran liblapack-dev + # glpk is not needed by Ipopt; it gives the certificate a second, + # independent opinion to be checked against. See + # TestAgreesWithLagrangian -- the Lagrangian spoke needs a solver that + # reports a dual bound, which Ipopt by definition does not, and which + # the bundle's cbc turns out not to do on this runner. + sudo apt-get install -y glpk-utils + + - name: Install Ipopt from idaes-ext + # The pip/conda ipopt is built against MUMPS only. The idaes-ext + # release bundle additionally carries the HSL linear solvers (ma27, + # ma57, ma97), which is what makes it usable on real NLPs. This is the + # same source Pyomo's and Egret's own CI pull from. + run: | + SOLVER_DIR="${GITHUB_WORKSPACE}/cache/solvers" + mkdir -p "$SOLVER_DIR" + echo "$SOLVER_DIR" >> $GITHUB_PATH + echo "LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$SOLVER_DIR" >> $GITHUB_ENV + URL=https://github.com/IDAES/idaes-ext + RELEASE=$(curl --max-time 150 --retry 8 -L -s \ + -H 'Accept: application/json' ${URL}/releases/latest) + VER=$(echo $RELEASE | sed -e 's/.*"tag_name":"\([^"]*\)".*/\1/') + # Fall back to a known-good release if the API is unavailable or + # rate-limited, so a GitHub hiccup does not read as a test failure. + if test -z "$VER" -o "$VER" = "$RELEASE"; then VER=3.4.2; fi + echo "idaes-ext version: $VER" + curl --max-time 300 --retry 8 -L \ + ${URL}/releases/download/$VER/idaes-solvers-ubuntu2204-x86_64.tar.gz \ + > solvers.tar.gz + tar -xzf solvers.tar.gz -C "$SOLVER_DIR" + "$SOLVER_DIR/ipopt" -v + + - name: setup the program + run: | + pip install -e ".[test]" + + - name: confirm the HSL linear solvers are really present + # Guards against silently falling back to a MUMPS-only build: a missing + # HSL solver makes ipopt exit abnormally rather than warn. + run: | + python -c " + import pyomo.environ as pyo + for ls in ('ma27', 'ma57', 'ma97'): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(-10, 10), initialize=0.0) + m.o = pyo.Objective(expr=(m.x - 3)**2) + m.c = pyo.Constraint(expr=m.x <= 1) + opt = pyo.SolverFactory('ipopt') + opt.options['linear_solver'] = ls + opt.solve(m) + assert abs(pyo.value(m.x) - 1.0) < 1e-6, ls + print(ls, 'OK') + " + echo "" + echo "==============================================================================" + echo "This build of Ipopt links HSL, and defaults to the ma27 linear solver, so" + echo "the tests below solve with HSL rather than with MUMPS." + echo "HSL, a collection of Fortran codes for large-scale scientific computation." + echo "See https://www.hsl.rl.ac.uk/" + echo "==============================================================================" + + - name: run ipopt-dependent tests + timeout-minutes: 10 + run: | + coverage run $COV_ARGS -m pytest -v \ + mpisppy/tests/test_dual_certificate.py \ + mpisppy/tests/test_ipopt_outer_bound.py + + - name: run the ipopt_outer_bound spoke on two ranks + timeout-minutes: 10 + run: | + mpiexec -np 2 coverage run $COV_ARGS -m mpi4py -m pytest -v \ + mpisppy/tests/test_ipopt_outer_bound.py + + - name: smoke-test the driver command line + # Stands in for a run_all.py entry: it exercises the same thing (the + # documented command line still works end to end) without forcing an + # ipopt install into the two run_all CI jobs, which have no other use + # for one. + timeout-minutes: 10 + run: | + cd examples/farmer + mpiexec -np 2 coverage run $COV_ARGS -m mpi4py \ + ../../mpisppy/generic_cylinders.py \ + --module-name farmer --num-scens 3 --solver-name ipopt \ + --max-iterations 5 --default-rho 1.0 --ipopt-outer-bound + + - name: List coverage files (debug) + if: always() + run: ls -la ${{ github.workspace }}/.coverage* || echo "No .coverage files found" + + - name: Upload coverage data + if: always() + uses: actions/upload-artifact@v4 + with: + name: coverage-ipopt-tests + path: ${{ github.workspace }}/.coverage.* + if-no-files-found: ignore + include-hidden-files: true + coverage-report: name: Coverage Report runs-on: ubuntu-latest @@ -1303,6 +1427,7 @@ jobs: - test-cg - test-agnostic - unit-tests + - ipopt-tests if: always() steps: - uses: actions/checkout@v3 diff --git a/.gitignore b/.gitignore index f5ea066e6..e609cddca 100644 --- a/.gitignore +++ b/.gitignore @@ -158,8 +158,18 @@ examples/**/*_full_solution/ examples/**/*_pickles/ examples/**/*_cyl_nonants.npy -# LaTeX build products under doc/slides (the .tex and the built .pdf are -# committed; everything else pdflatex/biber leaves behind is not) +# LaTeX build products under doc/slides and doc/designs (the .tex and the +# built .pdf are committed; everything else pdflatex/biber leaves behind is not) +doc/designs/**/*.aux +doc/designs/**/*.bbl +doc/designs/**/*.bcf +doc/designs/**/*.blg +doc/designs/**/*.fdb_latexmk +doc/designs/**/*.fls +doc/designs/**/*.out +doc/designs/**/*.run.xml +doc/designs/**/*.synctex.gz +doc/designs/**/*.toc doc/slides/**/*.aux doc/slides/**/*.bbl doc/slides/**/*.bcf diff --git a/doc/designs/ipopt_outer_bound_certificate.pdf b/doc/designs/ipopt_outer_bound_certificate.pdf new file mode 100644 index 000000000..327d4f293 Binary files /dev/null and b/doc/designs/ipopt_outer_bound_certificate.pdf differ diff --git a/doc/designs/ipopt_outer_bound_certificate.tex b/doc/designs/ipopt_outer_bound_certificate.tex new file mode 100644 index 000000000..91a1d45cc --- /dev/null +++ b/doc/designs/ipopt_outer_bound_certificate.tex @@ -0,0 +1,542 @@ +% Mathematical companion to doc/designs/ipopt_outer_bound_design.md +% Build: pdflatex ipopt_outer_bound_certificate.tex +\documentclass[11pt]{article} + +\usepackage[margin=1.1in]{geometry} +\usepackage{amsmath,amssymb,amsthm} +\usepackage{booktabs} +\usepackage[hidelinks]{hyperref} + +\theoremstyle{plain} +\newtheorem{lemma}{Lemma} +\newtheorem{theorem}[lemma]{Theorem} +\newtheorem{corollary}[lemma]{Corollary} +\newtheorem{proposition}[lemma]{Proposition} +\theoremstyle{remark} +\newtheorem{remark}[lemma]{Remark} + +\newcommand{\R}{\mathbb{R}} +\newcommand{\W}{W} +\newcommand{\vhat}{\hat{v}} +\newcommand{\qhat}{\hat{q}} +\newcommand{\OPT}{\mathrm{OPT}} + +\title{A certified outer bound for convex NLP scenario subproblems} +\author{mpi-sppy \\ \small companion to \texttt{doc/designs/ipopt\_outer\_bound\_design.md}} +\date{} + +\begin{document} +\maketitle + +\begin{abstract} +The usual device for getting outer bounds in mpi-sppy is to read the dual bound off the subproblem solver, which fails for Ipopt, which reports no dual bound. This note +derives the bound that the \texttt{ipopt\_outer\_bound} spoke computes instead. +The construction is Lagrangian weak duality combined with a tangent-plane +underestimator minimised in closed form over the variable box. \emph{None of the +mathematics is new}: Theorem~\ref{thm:main} is the Frank--Wolfe duality gap on a +box, and Section~\ref{sec:literature} places each result against the +literature. The derivation is written out because the hypotheses are what +matter in practice, and one of them is easy to get backwards. Its defining +property is that it requires \emph{no} assumption about the accuracy of the +solve: a truncated or sloppy solve yields a loose bound rather than an invalid +one, and an ill-conditioned one does the same, since the construction evaluates +rather than solves and so admits no amplification of error by a condition +number. Convexity, by contrast, is genuinely load-bearing, and +Section~\ref{sec:canonical} isolates the one place where it is easy to get +backwards. +\end{abstract} + +\section{Setting} + +We assume minimization. +Fix a scenario $s$ with probability $p_s$, $\sum_s p_s = 1$. Let +$v \in \R^{n}$ collect \emph{all} variables of the scenario subproblem --- the +nonanticipative variables $x(v)$ together with the recourse variables --- and let +\[ + B \;=\; \prod_{i=1}^{n} [\,l_i, u_i\,] +\] +be the box of variable bounds. With hub weights $\W_s$, the subproblem an outer bound +spoke solves is +\begin{equation}\label{eq:sub} + L_s(\W_s) \;=\; \min_{v} \;\bigl\{\, f_s(v) + \W_s^{\top} x(v) + \;:\; g_s(v) \le 0,\; h_s(v) = 0,\; v \in B \,\bigr\}, +\end{equation} +with $f_s : \R^n \to \R$, $g_s : \R^n \to \R^{m}$ and $h_s : \R^n \to \R^{k}$. +This is exactly the problem the ordinary Lagrangian spoke builds, with the +proximal term switched off. + +\paragraph{Why the returned objective value is not a bound here.} +For a minimisation, the value the solver reports is the objective evaluated at +the point it stopped at. Any such value is $\ge$ the minimum, so it is an +\emph{inner} bound: the wrong direction. Nor does Ipopt's convergence tolerance +close the gap, since it is a relative KKT-residual tolerance rather than a bound +on the objective error, and the two error sources point opposite ways --- +optimality error stops slightly above the optimum, while constraint violation +can place the iterate slightly below it. + +\section{Weak duality} + +\begin{lemma}[Weak duality]\label{lem:weak} +Let $\lambda \in \R^{m}_{+}$ and $\mu \in \R^{k}$ be arbitrary. Define +\begin{align} + \varphi_s(v) &\;=\; f_s(v) + \W_s^{\top} x(v) + + \lambda^{\top} g_s(v) + \mu^{\top} h_s(v), + \label{eq:phi}\\ + q_s(\lambda,\mu) &\;=\; \inf_{v \in B} \varphi_s(v). + \label{eq:q} +\end{align} +Then $q_s(\lambda,\mu) \le L_s(\W_s)$. +\end{lemma} + +\begin{proof} +Let $v$ be any point feasible for \eqref{eq:sub}. Then $v \in B$, so +$q_s(\lambda,\mu) \le \varphi_s(v)$. Moreover $g_s(v) \le 0$ and $\lambda \ge 0$ +give $\lambda^{\top} g_s(v) \le 0$, and $h_s(v) = 0$ gives +$\mu^{\top} h_s(v) = 0$; hence +$\varphi_s(v) \le f_s(v) + \W_s^{\top} x(v)$. Chaining the two inequalities, +$q_s(\lambda,\mu) \le f_s(v) + \W_s^{\top} x(v)$ for every feasible $v$. Taking +the infimum over feasible $v$ gives the claim. +\end{proof} + +Note what Lemma~\ref{lem:weak} does \emph{not} require: $\lambda$ and $\mu$ need +not be optimal, or even good. Any $\lambda \ge 0$ and any $\mu$ will do. This +is the property that makes the whole approach robust, and it is why a mistaken +sign convention or a stale multiplier can only cost tightness +(Remark~\ref{rem:robust}). + +\paragraph{The trap.} +Lemma~\ref{lem:weak} is not yet usable, because $q_s$ is defined by an +\emph{infimum}. Handing $\varphi_s$ to a nonlinear solver returns a +\emph{point}, and the value there is $\ge$ the infimum --- an upper bound on +$q_s$, which is the wrong direction again. A second NLP solve therefore does +not by itself certify anything. + +\section{The box underestimator} + +\begin{theorem}[Certified bound]\label{thm:main} +Suppose $\varphi_s$ is convex on an open set containing $B$ and differentiable +at $\vhat$. Define +\begin{equation}\label{eq:qhat} + \qhat_s(\vhat) \;=\; \varphi_s(\vhat) + \;+\; \sum_{i=1}^{n} \; \min_{t \in [l_i,u_i]} + \partial_i \varphi_s(\vhat)\,\bigl(t - \vhat_i\bigr). +\end{equation} +Then $\qhat_s(\vhat) \le q_s(\lambda,\mu) \le L_s(\W_s)$, and each term of the +sum is available in closed form: +\begin{equation}\label{eq:closed} + \min_{t \in [l_i,u_i]} \partial_i \varphi_s(\vhat)\,(t - \vhat_i) + \;=\; + \begin{cases} + \partial_i \varphi_s(\vhat)\,(l_i - \vhat_i), & \partial_i \varphi_s(\vhat) > 0,\\[2pt] + \partial_i \varphi_s(\vhat)\,(u_i - \vhat_i), & \partial_i \varphi_s(\vhat) < 0,\\[2pt] + 0, & \partial_i \varphi_s(\vhat) = 0. + \end{cases} +\end{equation} +\end{theorem} + +\begin{proof} +Convexity and differentiability at $\vhat$ give the gradient inequality +\[ + \varphi_s(v) \;\ge\; \varphi_s(\vhat) + + \nabla \varphi_s(\vhat)^{\top} (v - \vhat) + \qquad \text{for all } v . +\] +Minimising both sides over $v \in B$ preserves the inequality, and the right-hand +side is separable across coordinates because $B$ is a box, which yields +\eqref{eq:qhat}. Each one-dimensional problem minimises a linear function over +an interval, so the minimum is attained at an endpoint, giving +\eqref{eq:closed}. The final inequality is Lemma~\ref{lem:weak}. +\end{proof} + +\begin{corollary}\label{cor:anypoint} +$\qhat_s(\vhat) \le L_s(\W_s)$ for \emph{any} $\vhat$ at which $\varphi_s$ is +differentiable, \emph{any} $\lambda \ge 0$ and \emph{any} $\mu$. In particular +$\vhat$ need not be optimal, and need not even be feasible for \eqref{eq:sub}. +\end{corollary} + +Corollary~\ref{cor:anypoint} is the whole point. The certificate consumes the +iterate the solver happened to stop at and the multipliers it happened to +report, and returns a valid bound regardless. Computationally it costs one +gradient evaluation and a loop over the variables --- no second solve. It is +worth reading the corollary as covering \emph{ill-conditioned} solves and not +merely truncated ones; Remark~\ref{rem:conditioning} says why the distinction +does not arise here. + +\begin{remark}[Robustness]\label{rem:robust} +Because Corollary~\ref{cor:anypoint} holds for arbitrary $\lambda \ge 0$ and +$\mu$, an error in the sign convention used to recover the multipliers, a stale +multiplier, or a multiplier clipped to zero can only make $\qhat_s$ smaller. +Such an error makes the bound loose, never invalid. Convexity enjoys no such +protection; see Section~\ref{sec:canonical}. +\end{remark} + +\section{Exactness at a KKT point} + +The correction term in \eqref{eq:qhat} is not merely small in practice: it +vanishes identically at an exact KKT point, so the certificate degrades +gracefully and costs nothing when the solve is good. + +\begin{proposition}[The correction vanishes]\label{prop:kkt} +Let $\vhat$ satisfy the KKT conditions of \eqref{eq:sub} with multipliers +$\lambda \ge 0$ for $g_s \le 0$, $\mu$ for $h_s = 0$, and $z_L, z_U \ge 0$ for +the bounds $l - v \le 0$ and $v - u \le 0$. Then every term of the sum in +\eqref{eq:qhat} is zero, and if in addition $\varphi_s$ is convex then +$\qhat_s(\vhat) = L_s(\W_s)$. +\end{proposition} + +\begin{proof} +Stationarity of the full Lagrangian, which is $\varphi_s$ augmented with the +bound terms, reads +\begin{equation}\label{eq:stat} + \nabla \varphi_s(\vhat) \;=\; z_L - z_U . +\end{equation} +Fix a coordinate $i$ and consider the three possible positions of $\vhat_i$. + +If $l_i < \vhat_i < u_i$, complementarity gives $z_{L,i} = z_{U,i} = 0$, so +$\partial_i \varphi_s(\vhat) = 0$ by \eqref{eq:stat} and the term is $0$ by +\eqref{eq:closed}. + +If $\vhat_i = l_i$, complementarity gives $z_{U,i} = 0$, so +$\partial_i \varphi_s(\vhat) = z_{L,i} \ge 0$. When the derivative is strictly +positive \eqref{eq:closed} selects $t = l_i = \vhat_i$, and when it is zero the +term is zero outright; either way the term is $0$. + +If $\vhat_i = u_i$, symmetrically $z_{L,i} = 0$ and +$\partial_i \varphi_s(\vhat) = -z_{U,i} \le 0$, so \eqref{eq:closed} selects +$t = u_i = \vhat_i$ and the term is $0$. + +Hence $\qhat_s(\vhat) = \varphi_s(\vhat)$. Complementarity +$\lambda^{\top} g_s(\vhat) = 0$ and feasibility $h_s(\vhat) = 0$ reduce +\eqref{eq:phi} to $\varphi_s(\vhat) = f_s(\vhat) + \W_s^{\top} x(\vhat)$. Under +convexity the KKT conditions are sufficient for global optimality of +\eqref{eq:sub}, so that value is $L_s(\W_s)$. +\end{proof} + +Proposition~\ref{prop:kkt} has a practical corollary worth stating plainly: the +correction term \emph{measures its own inexactness}. It is zero when the solve +is exact and grows as the solve degrades, so the reported bound loosens smoothly +rather than becoming wrong. + +\section{Canonical form, and where convexity is easy to get backwards} +\label{sec:canonical} + +Theorem~\ref{thm:main} requires $\varphi_s$ convex, and by \eqref{eq:phi} that +requires each component of $g_s$ to be convex and each component of $h_s$ to be +affine, since $\lambda \ge 0$ but $\mu$ is free in sign. + +The subtlety is that $g_s$ is the constraint in \emph{canonical} form +$g_s(v) \le 0$, and putting a $\ge$ row into canonical form negates its body. +Writing $b(v)$ for the constraint body as the user stated it: + +\begin{center} +\begin{tabular}{lll} +\toprule +as written & canonical form & requirement on $b$ \\ +\midrule +$b(v) \le \beta$ & $g = b - \beta$ & $b$ convex \\ +$b(v) \ge \alpha$ & $g = \alpha - b$ & $b$ \textbf{concave} \\ +$\alpha \le b(v) \le \beta$ & both rows & $b$ affine \\ +$b(v) = \beta$ & $h = b - \beta$ & $b$ affine \\ +\bottomrule +\end{tabular} +\end{center} + +So $x^2 \le 4$ is admissible and $x^2 \ge 1$ is not, though both are written with +a convex body --- and indeed the feasible set of the second is not convex. This +is the one place in the construction where an entirely ordinary-looking model +silently leaves the hypotheses, and the consequence is not a loose bound but an +invalid one. A worked instance: for +\[ + \min\; x \quad \text{s.t.} \quad x^2 \ge 1, \quad x \in [0, 1.5], +\] +whose optimum is $1$, taking $\vhat = 0.5$ and $\lambda = 1$ in \eqref{eq:qhat} +yields $\qhat = 1.25 > 1$. + +Of the four rows above, three are decidable --- affineness is a question about +polynomial degree --- and the implementation enforces them. Convexity or +concavity of a one-sided nonlinear body is not decidable, and remains a user +assertion. + +\section{Aggregating across scenarios} + +The per-scenario bounds are combined by the usual probability-weighted sum, and +the condition under which that sum is itself a bound deserves care. + +\begin{proposition}[Aggregation]\label{prop:agg} +Let $\OPT$ be the optimal value of the original stochastic program. If the +weights satisfy $\sum_s p_s \W_s = 0$, then +$\sum_s p_s\, L_s(\W_s) \le \OPT$. +\end{proposition} + +\begin{proof} +Let $x^{*}$ together with the scenario recourse decisions be optimal for the +original problem, so that $\OPT = \sum_s p_s f_s(v_s^{*})$ where $v_s^{*}$ is the +optimal decision in scenario $s$ and $x(v_s^{*}) = x^{*}$ for every $s$ by +nonanticipativity. Each $v_s^{*}$ is feasible for subproblem \eqref{eq:sub}, so +$L_s(\W_s) \le f_s(v_s^{*}) + \W_s^{\top} x^{*}$. Multiplying by $p_s$ and +summing, +\[ + \sum_s p_s L_s(\W_s) + \;\le\; \sum_s p_s f_s(v_s^{*}) + \Bigl( \sum_s p_s \W_s \Bigr)^{\!\top} x^{*} + \;=\; \OPT + 0 . \qedhere +\] +\end{proof} + +\begin{remark}[Weights must share a generation]\label{rem:mixed} +The hypothesis $\sum_s p_s \W_s = 0$ is a \emph{joint} condition on the whole +weight vector, and progressive hedging maintains it at each iteration. It is +not preserved by mixing: if $\W'_s$ takes scenario $s$'s weights from one +iteration and scenario $t$'s from another, then in general +$\sum_s p_s \W'_s \ne 0$, the proof above leaves the residual term +$\bigl(\sum_s p_s \W'_s\bigr)^{\top} x^{*}$ in place, and the resulting number +is not a bound at all. Nothing downstream can detect this, which is why a +scenario whose certificate is unavailable must make the whole expectation +unavailable rather than contribute a stale value. +\end{remark} + +\section{Practical consequences}\label{sec:practical} + +\paragraph{Tightness is governed by the box.} +By \eqref{eq:closed} the shortfall of $\qhat_s$ below $\varphi_s(\vhat)$ is +\[ + \sum_{i=1}^{n} \bigl| \partial_i \varphi_s(\vhat) \bigr| \cdot + \bigl( \text{distance from } \vhat_i \text{ to the far end of } [l_i,u_i] \bigr), +\] +so the bound is only as tight as the variable bounds are. Gradient components +at a converged solve are small but generically nonzero, so a box of width $10$ +costs little while a box of width $10^{10}$ can cost a substantial fraction of +the objective. Tightening bounds before certifying --- for instance by +feasibility-based bounds tightening, which shrinks $B$ without removing feasible +points and therefore only raises $\qhat_s$ --- is the most effective way to +improve the bound. + +\begin{remark}[Admissible tightenings of the box]\label{rem:boxtighten} +The box in \eqref{eq:closed} need not be the one the model was written with, and +the implementation shrinks it by two different mechanisms. The condition under +which this is legitimate is not that a smaller box removes points --- that +argument shows the bound gets tighter, which is the direction that could break +it --- but the following. Let $x^{*}$ be an optimal solution of the full +problem. If $x^{*} \in B_s$ for every $s$, then for each $s$ +\[ + \qhat_s \;\le\; \inf_{v \in B_s} \varphi_s(v) \;\le\; \varphi_s(x^{*}) + \;\le\; f_s(x^{*}) + \W_s^{\top} x^{*}, +\] +the last step because $\lambda \ge 0$, $g_s(x^{*}) \le 0$ and $h_s(x^{*}) = 0$, +and Proposition~\ref{prop:agg} goes through unchanged. The requirement is +therefore that \emph{every $B_s$ contain one common optimal solution of the full +problem}; preserving a different optimiser in each scenario does not suffice, +since the weighted sum telescopes only at a single $x^{*}$. + +Two sufficient conditions are used. Feasibility-based bounds tightening +satisfies the requirement in the strong form: it removes no point satisfying the +constraints of subproblem $s$, so $B_s$ retains every feasible point of +\eqref{eq:sub} and in particular $x^{*}$, with no appeal to optimality. +Optimality-based tightening --- reduced-cost fixing, which does discard feasible +points --- satisfies only the weak form, and relies on its own contract that +some optimal solution survives, together with the fact that the same bounds are +applied to every scenario, so that the surviving solution is common to all of +them. A tightening that guarantees neither is not admissible, and no check in +the implementation can detect that. +\end{remark} + +\paragraph{Unbounded variables.} +If $l_i = -\infty$ and $\partial_i \varphi_s(\vhat) > 0$, or $u_i = +\infty$ and +$\partial_i \varphi_s(\vhat) < 0$, the corresponding term in \eqref{eq:closed} is +$-\infty$ and no finite bound is available. Whether this occurs depends on the +model and on how converged the solve is: at an exact KKT point the offending +component is zero by Proposition~\ref{prop:kkt}, but away from one it generally +is not. The implementation reports ``no bound'' in that case rather than +$-\infty$, so that Remark~\ref{rem:mixed} is respected. + +\paragraph{Inexact solves.} +Corollary~\ref{cor:anypoint} already covers these, but it is worth spelling out +what ``inexact'' can mean. Problem~\eqref{eq:sub} is convex by hypothesis, so +it has no non-global local minima. A solver returning a sub-optimal answer can +therefore only mean that it stopped short of converging, leaving an inexact +$\vhat$ and inexact multipliers. Both are admissible inputs to +Theorem~\ref{thm:main}: the conclusion $\qhat_s \le L_s(\W_s)$ is unchanged, and +the entire effect appears in the correction term of \eqref{eq:qhat}, which grows +as $\vhat$ moves away from a KKT point. By Proposition~\ref{prop:kkt} that term +vanishes exactly at a KKT point, so the correction is precisely the price of +inexactness, and it measures itself. + +\begin{remark}[Conditioning]\label{rem:conditioning} +The distinct worry is that \eqref{eq:sub} is badly conditioned, so that the +solve is not merely truncated but numerically poor. This affects the tightness +of $\qhat_s$ and not its validity, for a structural reason: \emph{the +certificate is an evaluation, not a solve.} A condition number quantifies how +much a linear solve amplifies error --- $\kappa(H)$ bounds the relative error in +$d$ obtained from $Hd = -g$ --- and no step of \eqref{eq:closed} inverts +anything. The computation evaluates $\varphi_s$ at $\vhat$, evaluates +$\nabla \varphi_s(\vhat)$, compares each component with zero, selects an +endpoint of the corresponding interval, multiplies and adds. Each is a forward +evaluation, and $\kappa$ has no route in. Correspondingly, the inequality the +construction rests on, +\[ + \varphi_s(v) \;\ge\; \varphi_s(\vhat) + + \nabla \varphi_s(\vhat)^{\top} (v - \vhat) + \qquad \text{for all } v, +\] +is a pointwise consequence of convexity. It holds exactly, at every $\vhat$, +with no error constant and no dependence on the conditioning of anything. +Conditioning governs how loose the plane is away from $\vhat$; it cannot change +which side of $\varphi_s$ the plane lies on. What it does change is the size of +$\nabla \varphi_s(\vhat)$ and the quality of the multipliers, and both push +$\qhat_s$ \emph{down}. +\end{remark} + +\paragraph{Floating point.} +Nothing in Theorem~\ref{thm:main} requires a tolerance. The arithmetic, being +carried out in finite precision, does introduce error, and the previous +paragraph says only that this error is not amplified by $\kappa$ --- not that it +is absent. + +\begin{remark}[Rounding]\label{rem:rounding} +Let $u$ denote the unit roundoff. Evaluating \eqref{eq:qhat} in floating point +incurs +\[ + \bigl| \mathrm{fl}(\qhat_s) - \qhat_s \bigr| + \;\lesssim\; + u \Bigl( |f_s(\vhat)| + + \textstyle\sum_i \lambda_i |g_i(\vhat)| + + \sum_j |\mu_j|\,|h_j(\vhat)| + + \sum_i T_i \, d_i \Bigr), +\] +where $T_i$ bounds the intermediate magnitudes arising in +$\partial_i \varphi_s(\vhat)$ and $d_i$ is the distance from $\vhat_i$ to the +far end of $[l_i, u_i]$. No condition number appears. +\end{remark} + +Remark~\ref{rem:rounding} is governed by the size of the \emph{summands}, so +cancellation matters: multipliers of order $10^{8}$ against an objective of +order one put the error floor near $10^{-8}$ even though $|\qhat_s|$ is of order +one. The implementation's cushion, reporting +$\qhat_s - \varepsilon(1 + |\qhat_s|)$, is scaled to the size of the +\emph{result} instead, and the two agree only on a well-scaled model. The +cushion is therefore a numerical safeguard and not part of the argument, and it +is exposed as a user-settable quantity rather than fixed. + +Note also that the final term of Remark~\ref{rem:rounding} is the one that can +raise $\qhat_s$. The correction in \eqref{eq:closed} equals +$-\sum_i |\partial_i \varphi_s(\vhat)|\, d_i$, so a gradient component computed +slightly small in magnitude makes it slightly less negative. Remark~\ref{rem:robust} +gives no protection here: it is a statement about the multipliers, and a +perturbation of $\nabla \varphi_s$ is not one-directional in the way a +perturbation of $\lambda$ is. + +\paragraph{Non-finite values.} +A solve that diverges can leave NaN in $\vhat$ or in the reported multipliers, +and an infinite multiplier yields NaN through $\infty \cdot 0$ or an infinite +$\varphi_s(\vhat)$. These propagate silently. The implementation therefore +rejects any non-finite $\qhat_s$ and reports ``no bound'', the same disposition +it takes for an unbounded direction, so that Remark~\ref{rem:mixed} is again +respected. The case that makes this necessary rather than cosmetic is +$+\infty$: a consumer keeping the best outer bound seen would discard NaN by +accident, since NaN fails every comparison, but would accept $+\infty$ as an +improvement. + +\section{Relation to known results}\label{sec:literature} + +Nothing in Sections~2--7 is new. This section attributes each result. + +\paragraph{Lemma~\ref{lem:weak}} is textbook Lagrangian weak duality, stated for +an arbitrary multiplier pair rather than an optimal one. + +\paragraph{Theorem~\ref{thm:main}} is the \emph{Frank--Wolfe duality gap}, +sometimes called the Wolfe gap or the linearisation gap. For a convex $f$ on a +compact convex $D$, evaluating the linear minimisation oracle at $x$ gives +\[ + f(x) + \min_{y \in D} \nabla f(x)^{\top}(y - x) \;\le\; \min_{D} f , +\] +which is \eqref{eq:qhat} with $D = B$. The bound is used in the Frank--Wolfe +algorithm \cite{frankwolfe1956}, where the same oracle evaluation that produces +it also supplies the search direction; reading it instead as a +\emph{certificate} computable at any iterate is the framing in +\cite{jaggi2013}. The only specialisation here is that a box makes the oracle +separable and closed-form, so the certificate costs one gradient evaluation and +no optimisation at all. + +\paragraph{Proposition~\ref{prop:agg}} is the standard argument by which +progressive hedging weights yield a lower bound on the optimal value, as in +\cite{gade2016}. + +\paragraph{The remarks of Section~\ref{sec:practical}} are elementary rather +than novel. Remark~\ref{rem:rounding} is the textbook forward-error bound for a +floating-point sum, specialised to the summands of \eqref{eq:qhat}, and +Remark~\ref{rem:boxtighten} is two applications of Lemma~\ref{lem:weak}. They +are stated because the hypotheses they turn on are easy to lose sight of in +implementation, not because they are results. + +\paragraph{The combination} of progressive hedging with Frank--Wolfe machinery to +obtain Lagrangian dual bounds is itself established \cite{boland2018}, and +mpi-sppy already ships an FWPH cylinder on that basis. + +What is specific here is therefore an assembly rather than a result: the +multipliers are the ones Ipopt already returns, the linearisation point is the +iterate it already stopped at, and the box is the model's own variable bounds --- +so a valid outer bound is obtained from the solve the spoke was going to perform +anyway, with no second optimisation and no assumption that the first one +converged. The wider question of extracting valid dual bounds from an inexact +oracle is treated at length in the inexact bundle-method literature, which this +note does not attempt to survey. + +\section{Correspondence with the implementation} + +\begin{center} +\begin{tabular}{ll} +\toprule +object here & in the code \\ +\midrule +$\varphi_s$ in \eqref{eq:phi} & \texttt{\_lagrangian\_expression} \\ +$\qhat_s$ in \eqref{eq:qhat} & \texttt{certified\_lower\_bound} \\ +decidable hypotheses & \texttt{check\_model\_is\_certifiable} \\ +box tightening, strong form & \texttt{unbounded\_variables} (\texttt{fbbt}) \\ +Remark~\ref{rem:rounding} cushion & \texttt{eps\_rel} \\ +non-finite rejection & \texttt{math.isfinite} test on $\qhat_s$ \\ +\midrule +box tightening, weak form & \texttt{receive\_nonant\_bounds} \\ +Proposition~\ref{prop:agg} sum & \texttt{SPOpt.Ebound} \\ +Remark~\ref{rem:mixed} safeguard & a \texttt{None} bound, propagated collectively \\ +\bottomrule +\end{tabular} +\end{center} + +\noindent +Everything above the rule is in \texttt{utils/dual\_certificate.py}, +whose \texttt{eps\_rel} argument is exposed on the command line as +\texttt{--ipopt-outer-bound-cushion}. Below it, the weak-form tightening is in +\texttt{cylinders/spcommunicator.py} and the sum is in \texttt{spopt.py}, with +the \texttt{None} propagated across both. All paths are relative to +\texttt{mpisppy/}, and the spoke driving the whole is +\texttt{cylinders/ipopt\_outer\_bound.py}. + +\begin{thebibliography}{9} + +\bibitem{frankwolfe1956} +Frank, M. and Wolfe, P.: +\newblock An algorithm for quadratic programming. +\newblock \emph{Naval Research Logistics Quarterly} \textbf{3}(1--2), 95--110 (1956). + +\bibitem{jaggi2013} +Jaggi, M.: +\newblock Revisiting Frank--Wolfe: projection-free sparse convex optimization. +\newblock \emph{Proceedings of the 30th International Conference on Machine +Learning (ICML)}, 427--435 (2013). + +\bibitem{gade2016} +Gade, D., Hackebeil, G., Ryan, S.M., Watson, J.-P., Wets, R.J.-B. and +Woodruff, D.L.: +\newblock Obtaining lower bounds from the progressive hedging algorithm for +stochastic mixed-integer programs. +\newblock \emph{Mathematical Programming} \textbf{157}(1), 47--67 (2016). + +\bibitem{boland2018} +Boland, N., Christiansen, J., Dandurand, B., Eberhard, A., Linderoth, J., +Luedtke, J. and Oliveira, F.: +\newblock Combining progressive hedging with a Frank--Wolfe method to compute +Lagrangian dual bounds in stochastic mixed-integer programming. +\newblock \emph{SIAM Journal on Optimization} \textbf{28}(2), 1312--1336 (2018). + +\end{thebibliography} + +\end{document} diff --git a/doc/designs/ipopt_outer_bound_design.md b/doc/designs/ipopt_outer_bound_design.md new file mode 100644 index 000000000..45057c7b6 --- /dev/null +++ b/doc/designs/ipopt_outer_bound_design.md @@ -0,0 +1,806 @@ +# `ipopt_outer_bound` — a certified outer-bound cylinder for convex NLP subproblems + +Status: draft for review. Branch `ipopt-outer-bound` (off Pyomo/mpi-sppy `main`). + +The mathematics is set out separately, with proofs, in +[`ipopt_outer_bound_certificate.tex`](ipopt_outer_bound_certificate.pdf) — +weak duality, the box underestimator, exactness at a KKT point, the canonical-form +sign condition, and the aggregation hypothesis. This document covers the design +decisions and the implementation; that one covers why the bound is valid. + +## 1. Goal + +Give mpi-sppy an outer bound for stochastic programs whose scenario subproblems are +**convex NLPs solved with Ipopt**. Today there is none: the Lagrangian spoke's bound +comes from the solver's dual bound, Ipopt is not a branch-and-bound solver and reports +no dual bound, and the spoke says so out loud: + +```python +# mpisppy/cylinders/lagrangian_bounder.py +if "ipopt" in self.opt.options["solver_name"]: + print("\n WARNING: An ipopt solver will not give outer bounds\n") +``` + +Measured on this branch, a converged Ipopt solve through Pyomo returns: + +| field | value | +|---|---| +| `results.problem.lower_bound` / `upper_bound` | `-inf` / `inf` (Pyomo's untouched defaults) | +| `results.solution[0].objective` | `{}` — empty, even with `load_solutions=False` | +| `results.solver.message` | `'Ipopt 3.13.2\x3a Optimal Solution Found'` | + +So `Ebound()` sums `-inf` and a user running a convex stochastic NLP gets no bound at +all. This design supplies one that is *certified* — valid by a theorem, not by an +assertion that the solver converged. + +That phrase is meant literally, and the obvious objection to it — *Ipopt is an +iterative method on a possibly ill-conditioned problem, so how can its output certify +anything?* — is answered in **§5.3**. The short form: the certificate is an evaluation, +not a solve, so there is no linear system for a condition number to amplify, and the +inequality it rests on holds exactly at every point. Ill-conditioning costs tightness. +The two places where numerics do bear on validity, both of them about scaling rather +than conditioning, are in §5.2. + +Scope is **Ipopt only**, by decision. The mechanism generalizes to any NLP solver +returning duals, but nothing here is written to be solver-neutral, and the dual sign +conventions in §5 are measured from Ipopt. + +## 2. Why this cannot just read a number off the solver + +For a minimization subproblem the solver's returned objective value is an *inner* +bound: it is the value at a point, hence ≥ the subproblem optimum. An outer bound +needs a number ≤ the optimum. Ipopt supplies the former and never the latter. + +Ipopt's `tol` does not close the gap either. It is a **relative KKT-residual** +tolerance (the binary's own `-=` listing: "Desired convergence tolerance (relative)"), +not an objective-error tolerance; converting one to the other needs multiplier +magnitudes and curvature. And the two error sources point opposite ways: optimality +error stops you slightly *above* the optimum (unsafe for an outer bound), while +constraint violation — `constr_viol_tol` defaults to `1e-4`, four orders looser than +`tol` — can put you slightly *below* it. A single scalar cushion covers neither +honestly. + +What *is* rigorous is Lagrangian weak duality, which needs no convergence assumption +at all. That is what this cylinder computes. + +## 3. The bound + +### 3.1 Setting + +The hub relaxes non-anticipativity with weights `W`. Provided the weights satisfy the +usual condition `Σ_s p_s W_s = 0`, the scenario-separable Lagrangian dual value + +``` + D(W) = Σ_s p_s L_s(W_s) ≤ optimum +``` + +is an outer bound, where the scenario subproblem — exactly the problem the Lagrangian +spoke already builds and solves, prox term off — is + +``` + L_s(W_s) = min f_s(v) + W_sᵀ x(v) + s.t. g_s(v) ≤ 0, h_s(v) = 0, v ∈ B = [lo, hi] +``` + +`v` is every variable of the scenario model (non-anticipative `x` plus recourse), and +`B` is its box of variable bounds. + +**Convexity assumption (load-bearing, unverifiable in general):** `f_s` convex, each +component of `g_s` convex, `h_s` affine, all over `B`. §6 lists the parts of this that +*are* mechanically checkable; the rest is a user assertion. + +Note that the requirement is on the **canonical** `g`, which negates the body of a `>=` +row, so it is not the same as "the constraint body is convex": + +| as written | canonical `g` | requirement on the body | +|---|---|---| +| `body ≤ upper` | `body − upper` | convex | +| `body ≥ lower` | `lower − body` | **concave** | +| `lo ≤ body ≤ up` | both rows | affine | +| `body == rhs` | `body − rhs` | affine | + +The theorem applies to `x² ≤ 4` but not to `x² ≥ 1`, though both are written with a +convex body — and the feasible set of the second is not convex at all. This is worth +stating loudly because it is easy to read the wrong way: a code review of this branch +turned up that the guard, the unit test and `spokes.rst` had all settled on "only +equalities need to be affine", which let `min x s.t. x² ≥ 1, x ∈ [0, 1.5]` through and +certified 1.25 for a problem whose optimum is 1.0 — an outer bound above the optimum. + +### 3.2 The dual function and the trap + +For multipliers `λ ≥ 0` and free `μ`, define + +``` + φ_s(v) = f_s(v) + W_sᵀ x(v) + λᵀ g_s(v) + μᵀ h_s(v) + q_s(λ, μ) = inf_{v ∈ B} φ_s(v) ≤ L_s(W_s) (weak duality) +``` + +Weak duality holds for **any** `λ ≥ 0` and **any** `μ` — no optimality of the +multipliers is required. That is the whole reason this approach is rigorous where +"trust the objective value" is not. + +The trap: `q_s` is defined by an *infimum*. Handing `φ_s` to Ipopt and solving it +returns a point, and the value at that point is ≥ the infimum — an upper bound on +`q_s`, which is the wrong direction again. **A second NLP solve does not by itself +produce a certificate.** This is the one place where the original sketch for this +cylinder (solve, take λ, solve the dual, report) does not close. + +### 3.3 What does close it: a linear underestimator over the box + +`φ_s` is convex on `B`, so for any point `v̂ ∈ B` the tangent at `v̂` lies below it: + +``` + φ_s(v) ≥ φ_s(v̂) + ∇φ_s(v̂)ᵀ (v − v̂) for all v +``` + +Minimizing the right-hand side over the box is separable and closed-form, so + +``` + q_s(λ, μ) ≥ φ_s(v̂) + Σ_i min_{v_i ∈ [lo_i, hi_i]} ∂_i φ_s(v̂) · (v_i − v̂_i) + ╰──────────────────────────────────────────────────╯ + = ∂_i φ · (lo_i − v̂_i) if ∂_i φ > 0 + ∂_i φ · (hi_i − v̂_i) otherwise +``` + +Call the right-hand side `q̂_s`. Then + +``` + q̂_s ≤ L_s(W_s) for ANY v̂ ∈ B, ANY λ ≥ 0, ANY μ +``` + +and the cylinder reports `Σ_s p_s q̂_s`. One gradient evaluation and a loop over +variables — no second solve, no tolerance argument, no feasibility requirement on `v̂`. + +**Looseness has a closed form, and it is the box width that sets it.** Each term is +`|∂_i φ|` times the distance from `v̂_i` to the far end of its interval, so + +``` + L_s(W_s) − q̂_s ≈ Σ_i |∂_i φ(v̂)| · (width of the box in the descending direction) +``` + +Confirmed numerically: on the running example a converged solve leaves +`|∂_x φ| = 5.1e−11` over a box of width 9, and the observed gap to the analytic optimum +is 4.56e−10 ≈ 5.1e−11 × 9. The practical consequence is the one worth telling users: +**this cylinder is only as tight as the variable bounds are.** A model whose bounds come +back from `fbbt` as ±1e10 will produce valid but useless numbers — 5e−11 × 1e10 is half a +unit of objective — which is the same failure the unbounded case in §6.1 reaches in the +limit, not a different one. + +### 3.4 Why the correction term is ~0 in normal operation + +At an exact KKT point of the scenario subproblem, with its multipliers `(λ, μ)` and +bound multipliers `z_L, z_U ≥ 0`, stationarity says exactly + +``` + ∇φ_s(v̂) = z_L − z_U +``` + +Componentwise: an interior `v̂_i` has `z_L,i = z_U,i = 0`, so `∂_i φ = 0` and the term +is 0. A `v̂_i` sitting at `lo_i` has `∂_i φ = z_L,i ≥ 0`, so the minimizing `v_i` is +`lo_i = v̂_i` and the term is again 0; symmetrically at `hi_i`. **The correction +vanishes at an exact KKT point**, where `q̂_s = φ_s(v̂) = f_s(v̂) + W_sᵀx̂ = L_s(W_s)` +by complementarity — strong duality, recovered exactly. + +So the correction term is precisely the price of inexactness, and it measures itself. +Measured on a 3-variable convex NLP with one convex `≤`, one active `≥`, one equality, +and finite bounds: + +| solve | `φ(v̂)` | correction | certified bound | vs. reference | +|---|---|---|---|---| +| converged | 15.93325207 | −3.6e−09 | 15.93325207 | tight to 9 digits | +| `max_iter=5` | 15.93321 | −3.8e−02 | 15.89485 | valid, 0.04 loose | +| `max_iter=3` | 14.35348 | −2.6e+01 | −11.85 | valid, useless | +| `max_iter=2` | 18.86760 | −6.9e+01 | −49.97 | valid, useless | + +Note the `max_iter=2` row: the naive "use the objective value" bound would have been +18.87, which is **above** the optimum — an invalid outer bound. The certificate +returns −49.97 instead: worthless, but sound. That is the trade this design makes +everywhere. + +### 3.5 Consequence for cost: one solve, not two + +Because `v̂` from the ordinary subproblem solve is already the minimizer of `φ_s` over +the box (§3.4), the anticipated second solve per iteration is very nearly a no-op at a +converged first solve. **The cylinder costs the same as the existing Lagrangian spoke +— one solve per scenario per iteration** — and lags the hub no more than that spoke +does. The re-solve retains value only when the first solve is sloppy, so it becomes an +opt-in tightening pass (Phase 4) rather than the mechanism. + +### 3.6 Which tightenings of the box are admissible + +The box is not fixed. Two separate mechanisms shrink it, they shrink it for different +reasons, and they are **not** sound for the same reason — which is worth stating, +because the tempting one-line justification ("tightening the box only removes points, +and fewer points can only raise an infimum") is an argument that the bound gets +*tighter*, not an argument that it stays *valid*. Raising `q̂_s` is exactly the +direction that could break it. + +What the aggregate bound actually needs is easy to state. Let `x*` be an optimal +nonanticipative solution of the full problem. For each scenario `s`, if `x* ∈ B_s` +then + +``` + q̂_s ≤ inf_{v ∈ B_s} φ_s(v) ≤ φ_s(x*) ≤ f_s(x*) + W_sᵀx* +``` + +— the last step because `λ ≥ 0`, `g_s(x*) ≤ 0` and `h_s(x*) = 0` — and summing with +weights `p_s` under `Σ_s p_s W_s = 0` gives `Σ_s p_s q̂_s ≤ OPT`, which is §8. So the +requirement on the boxes is precisely: + +> **every `B_s` must contain one common optimal solution of the full problem.** + +Note "common": preserving a *different* optimizer in each scenario would not do, since +the sum telescopes only at a single `x*`. + +**fbbt (setup, once).** `unbounded_variables(s, do_fbbt=True)` runs feasibility-based +bounds tightening before any certificate is computed. This satisfies the requirement +in its strong form: fbbt removes no point that satisfies the scenario's constraints, +so `B_s` still contains *every* feasible point of subproblem `s`, `x*` among them. No +appeal to optimality is needed, which is why this one is safe to describe as "shrinks +`B` without removing a feasible point". + +**The nonant-bounds channel (every iteration).** `receive_nonant_bounds()` narrows the +nonanticipative variables' bounds from `Field.NONANT_LOWER/UPPER_BOUNDS`, and in the +current code base only `reduced_costs_spoke` sends that field. Reduced-cost fixing is +*optimality*-based: it discards points that are provably no better than an incumbent, +which does remove feasible points. It therefore does **not** satisfy the strong form, +and the weak form is what carries it — its contract is that at least one optimal +solution survives, and because the bounds are broadcast and applied identically to +every scenario's nonants, the solution that survives is the same one everywhere. That +is the "common `x*`" the requirement asks for. + +Two practical notes. First, this is close to moot today: reduced-cost fixing wants +discrete variables, which this cylinder rejects as a hard error at setup, so in +practice the field is not being sent to it. Second, fbbt is deliberately not re-run +after nonant bounds arrive. That leaves tightness on the table rather than risking +anything — re-running it would still be sound, since anything fbbt infers from +constraints and bounds that all hold at `x*` also holds at `x*` — but the setup-time +pass is where the unbounded-variable warning in §6.1 belongs, and repeating it every +iteration would buy little. + +**Anything else that narrows the box needs one of these two arguments made for it.** +A tightening that preserves neither all feasible points nor a common optimum breaks the +bound, and nothing in the code can detect that. + +## 4. Getting the gradient + +`∇φ_s(v̂)` comes from Pyomo's reverse-mode differentiation: + +```python +from pyomo.core.expr.calculus.derivatives import differentiate, Modes +grad = differentiate(phi, wrt_list=vlist, mode=Modes.reverse_numeric) +``` + +Measured on this branch, a 500-variable / 500-constraint model: **7.6 ms** for the +whole gradient — negligible against an NLP solve, and one reverse sweep rather than +one pass per variable. + +This deliberately avoids PyNumero. `PyomoNLP` would also serve, but it needs the +compiled `pynumero_ASL` library, and `AmplInterface.available()` is `False` here — it is +not an mpi-sppy dependency and nothing installs it. The `differentiate` route adds no +dependency at all. + +Worth recording precisely, because it is easy to misread: the idaes-ext bundle §10.1 +installs *does* ship `libpynumero_ASL.so`. That does not make PyNumero available, because +Pyomo looks for it under `PYOMO_CONFIG_DIR` (`~/.pyomo`), not on `PATH`. So the choice +here is not "PyNumero is unobtainable" — it is one deliberate copy away — but rather that +taking it would add an install step and a second way for the gradient to be unavailable +at runtime, in exchange for nothing this calculation needs. + +## 5. Multipliers: canonical form and Ipopt's sign conventions + +Constraints are canonicalized so every inequality reads `g(v) ≤ 0`. The mapping from +Pyomo's `dual` suffix to the canonical multiplier is orientation-dependent; measured +against analytically-known multipliers (`min (x−3)²` with the constraint active, true +multiplier 4): + +| Pyomo constraint | canonical form | Pyomo `dual` | canonical multiplier | +|---|---|---|---| +| `body <= upper` | `g = body − upper` | −4.0 | `λ = max(−d, 0)` | +| `body >= lower` | `g = lower − body` | +4.0 | `λ = max(+d, 0)` | +| `body == rhs` | `h = body − rhs` | −4.0 | `μ = −d` | + +The equality row agrees with the existing repo constant +`solver_dual_sign_convention['ipopt'] = -1` in `mpisppy/utils/lshaped_cuts.py`. + +Ranged constraints (`lower ≤ body ≤ upper`, both finite) carry one dual for two rows; +they are split into both inequalities and the two rules above are applied, which lands +the magnitude on the active side and zero on the other. Verified in both directions: a +range active at its upper bound returns `d = −4` (`λ_upper = 4`, `λ_lower = 0`), and one +active at its lower bound returns `d = +4` (`λ_lower = 4`, `λ_upper = 0`). An inactive +inequality returns `d = 0`, so both clipped multipliers are 0 and the constraint drops +out of `φ` — which is what complementarity requires. + +Bound multipliers `z_L`/`z_U` are **not needed** — the box is kept explicit in the +certificate rather than dualized, so `ipopt_zL_out`/`ipopt_zU_out` need not be +imported at all. + +### 5.1 The robustness property that makes this safe + +Weak duality holds for any `λ ≥ 0` and any `μ`. Therefore **a sign-convention error, a +stale dual, or a clipped multiplier can only make the bound loose, never wrong.** The +observed failure mode of a wrong sign is a bound of −34.7 where 1.2476 was available: +obviously broken, harmless. Combined with the hub keeping the best outer bound seen +(`_outer_bound_update` in `spcommunicator.py`), a weak bound is silently ignored. + +The exception is convexity, which is genuinely load-bearing: if the model is not +convex, the tangent in §3.3 is not an underestimator and the bound is simply wrong. +Hence the guards in §6. + +### 5.2 Floating point + +At a converged solve the certified bound exceeded `f(v̂)` by 1.5e−06. This is not a +demonstrated violation — `v̂` carries up to `constr_viol_tol` of infeasibility, so +`f(v̂)` understates the true optimum — but it shows the margin sits at the level of the +solver's feasibility tolerance. A small relative cushion is cheap insurance and costs +nothing that matters in a PH outer bound. + +**Decided: the cushion is on by default.** The engine reports + +``` + q̂ − ε_rel·(1 + |q̂|), ε_rel = 1e-9 by default +``` + +settable through `--ipopt-outer-bound-cushion` (0 disables it). Note honestly what this +does and does not buy: 1e-9 is an order of magnitude *smaller* than the 1.5e−06 margin +measured above, so it is last-bit hygiene, not a proof-carrying margin. The margin +itself is not a soundness problem — the theorem in §3.3 needs no tolerance argument — +and a user who wants to absorb `constr_viol_tol`-scale infeasibility can raise `ε_rel`. + +**Where the cushion does not scale with the risk.** `φ(v̂)` is evaluated in double +precision as `f + Σᵢ λᵢgᵢ + Σⱼ μⱼhⱼ`, and `q̂` adds `−Σᵢ |∂ᵢφ|·dᵢ`, where `dᵢ` is the +distance from `v̂ᵢ` to the far end of its interval. Rounding gives + +``` + |q̂_computed − q̂_exact| ≲ u·( |f| + Σᵢ|λᵢgᵢ| + Σⱼ|μⱼhⱼ| ) + u·Σᵢ (terms of ∂ᵢφ)·dᵢ + u = 2⁻⁵³ ≈ 1.1e−16 +``` + +Both error terms are governed by the size of the **summands**. The cushion +`ε_rel·(1+|q̂|)` is governed by the size of the **result**. On a well-scaled model those +track each other and 1e−9 is generous by seven orders of magnitude. Where cancellation +is severe they come apart: a constraint row scaled by 1e−8 drives its multiplier to +~1e8 against an objective of order one, so the error floor is ~1e−8 absolute while the +cushion is ~1e−9. + +Note that the second error term is the one that can push `q̂` *up*. The correction is +`−Σᵢ|∂ᵢφ|dᵢ`, so a gradient component computed slightly too small in magnitude makes +the correction slightly less negative. Unlike a multiplier error, which §5.1 shows is +one-directional, a gradient error is not. It is not amplified by anything — see §5.3 — +but it is not self-protecting either. + +**Practical consequence: on a model known to be badly scaled, raise +`--ipopt-outer-bound-cushion`.** This is the one respect in which numerics bear on +validity rather than on tightness, and it is the reason the flag exists as a knob +instead of a constant. + +**Non-finite results are screened.** A diverged solve can leave NaN in the point or in +the duals, and NaN propagates silently through every expression above. An infinite +multiplier produces NaN (`inf·0`) or `±inf`. `certified_lower_bound` therefore checks +`math.isfinite(q̂)` before applying the cushion and returns `None` otherwise, reusing the +word it already has for "no bound this time". Without the check the safety would be +accidental: NaN loses every comparison, so the hub's `new > old` test in +`_outer_bound_update` happens to reject it — but `+inf` would compare as an +*improvement* and be latched as an outer bound that is not one. + +### 5.3 Ill-conditioning: why it costs tightness and not validity + +Ill-conditioning is the obvious worry for a bound read out of an NLP solver, so it is +worth saying precisely where it does and does not enter. + +**"Sub-optimal" has only one meaning here.** The model is convex by hypothesis, so it +has no non-global local minima. Ipopt returning a sub-optimal answer can therefore only +mean it stopped short of converging — an inexact iterate with inexact multipliers. That +is exactly the case §3.3 was constructed to cover: the corollary holds for *any* `v̂`, +optimal or not, feasible or not, and for *any* `λ ≥ 0` and *any* `μ`. + +**The bounding plane is immune because there is no solve.** A condition number measures +how much a *linear solve* amplifies error — `κ(H)` bounds the relative error in `d` +from `Hd = −g`. The certificate inverts nothing. Its entire computation is: evaluate +`φ` at `v̂`, evaluate `∇φ` at `v̂` by reverse-mode AD, compare each component to zero, +select an endpoint of that variable's interval, multiply, and add. Every one of those +is a forward evaluation carrying relative error `O(u)` per operation. `κ` has no route +in. Correspondingly, the inequality the plane rests on, + +``` + φ(v) ≥ φ(v̂) + ∇φ(v̂)ᵀ(v − v̂) for all v, +``` + +is a pointwise consequence of convexity. It holds exactly, at every `v̂`, with no error +constant and no dependence on how well-conditioned anything is. Conditioning changes +how *loose* the plane is away from `v̂`; it cannot change which side of `φ` it lies on. + +**Where conditioning does land: the correction term.** A badly conditioned problem +leaves `v̂` further from optimal and `∇φ(v̂)` correspondingly larger, and the shortfall +`Σᵢ|∂ᵢφ|·dᵢ` grows with it. It also gives worse multipliers, which weakens `φ` itself. +Both effects move the bound *down*. This is visible rather than theoretical: in the +study below the well-scaled Hilbert QP closed to 1e−9 relative, while the same problem +with a 1e−8 row scaling stalled at 1.8e−2 relative and never closed further. Loose, +valid, and correctly reported as such. + +**Empirical check.** Hilbert-matrix QPs (`H_ij = 1/(i+j+1)`, `cond(H)` ≈ 1.6e13 at +n=10 and 3.5e17 at n=16), with row scalings of 1e±8 and an objective offset of 1e12, +solved at `max_iter ∈ {1,2,3,5,10,∞}` and certified at each. Every bound was compared +against `f` at a point *constructed* to be feasible — not against `f(v̂)`, which is the +trap: `v̂` carries up to `constr_viol_tol` of infeasibility, so `f(v̂)` can sit below the +true optimum and manufacture an apparent violation where there is none. No violations at +any conditioning, scaling, or truncation level. This is `TestIllConditioning` in +`mpisppy/tests/test_dual_certificate.py`. + +**What ill-conditioning *can* break is convexity, not the certificate.** A Hessian that +is nearly singular and slightly indefinite is non-convex, the tangent is then not an +underestimator, and §5.1's exception applies with full force. That is a property of the +model as written, and no guard in §6 catches it. + +## 6. Guards + +Checkable at setup, hard error (following the repo's fail-loudly convention): + +- **Any integer or binary variable** in a subproblem. A convexity claim is + definitionally false there, and this is the failure mode most likely to be reached by + accident. +- **Any equality constraint with `polynomial_degree() != 1`.** A nonlinear equality + makes `μᵀh` non-convex for one sign of `μ`, breaking §3.3. +- **Any two-sided (ranged) constraint with `polynomial_degree() != 1`.** It splits into + `g = body − upper` *and* `g = lower − body`, so its body would have to be both convex + and concave — affine. Like the equality case this is decidable, so it is enforced + rather than assumed. One-sided nonlinear rows are *not* rejected: whether the body is + convex (for `≤`) or concave (for `≥`) is not decidable here, and is the user's + assertion. See the table in §3.1 — the direction of that assertion flips with the + orientation, which is the part most likely to be got wrong. +- **Prox term attached.** With a prox term the subproblem is not the Lagrangian + relaxation and the bound is not a Lagrangian bound. `lagrangian_prep` already passes + `attach_prox=False`; assert it. +- **Maximization.** Minimize-only, with a hard error. The concave mirror that would add + max support is deferred and may not happen; see §10. +- **Solver is not Ipopt.** Scope decision; the sign conventions in §5 are measured from + Ipopt only. +- **A `dual` Suffix that does not import.** A scenario creator may already have attached + one — supplying dual warm starts is a common reason — and silently reusing an + `EXPORT` or `LOCAL` suffix would import nothing, leaving the certificate with no + multipliers. Existence is therefore not enough; `import_enabled()` is what is checked. + +Checked every iteration rather than at setup, and **not** a hard error — warn once on +rank 0 and report no bound, for the reason in §6.1: + +- **Nonanticipative variables fixed after setup.** Same failure as the prox term and + less obvious: PH's variable-fixing extensions (`fixer.py`, the reduced-cost fixers) + restrict the subproblem, which can only *raise* its minimum, so the resulting number + bounds the restricted problem and not the original. The certificate engine cannot see + this — a fixed variable is simply a constant to it — so the cylinder snapshots which + nonants were fixed at setup and compares each iteration. Nonants the *scenario + creator* fixed are part of the problem and are fine; only newly fixed ones suppress + the bound. This is a real interaction: fixing extensions are commonly on. It is not a + hard error because the fixing may be transient, and because raising from inside the + iteration loop would `MPI_Abort` the hub and every other spoke. + +### 6.1 Unbounded variables: warn and stand down, do not fail + +The closed-form box minimization in §3.3 returns −∞ as soon as an unbounded direction +has a nonzero gradient component. + +How often that happens is **model-dependent, and the first draft of this design +overstated it.** Measured on the running example at a converged solve, one component is +−5.1e−11 (nonzero, as claimed) while the component belonging to the variable that was +left unbounded is *exactly* 0.0 — the cancellation `2(ŷ−2) + λ` is exact there — so a +bound is still available despite the missing bounds. Truncate the same solve to one +iteration and that component becomes −1.9e−01 and the certificate correctly returns +`None`. The honest statement: at a converged solve these components are small but +generically nonzero, so an unbounded variable *may* defeat the certificate and reliably +does whenever the solve is not converged. Neither outcome can be counted on in advance, +which is what makes warn-and-stand-down the right policy rather than either +fail-at-setup or assume-it-is-fine. + +`fbbt` runs first (precedent: the CVaR `eta` bound in `utils/cvar.py`), which recovers +bounds implied by the constraints. **Decided: a variable still unbounded after `fbbt` +does not raise.** Setup emits a warning naming the offending variables, and the +cylinder then reports no bound — `outer_bound = None` per scenario, so `Ebound()` +returns `None` collectively and nothing is sent to the hub. + +Rationale for warning rather than erroring: this cylinder is an *optional* source of a +bound. A model that is merely under-bounded is not a broken model, and killing a whole +parallel run over a spoke that could have sat quiet is the wrong trade — unlike the +guards above, where the model genuinely violates an assumption and any number the +cylinder produced would be wrong. The warning is what keeps the quiet case from being +silent. + +Per repo convention the warning is emitted on one rank only (`cylinder_rank == 0`), +once at setup rather than once per iteration. + +Not checkable, documented as the user's assertion: convexity of `f_s` and `g_s`. + +## 7. Solver scoping: Ipopt for this cylinder, anything for everyone else + +This cylinder pins **its own** solver to Ipopt. The hub and every other spoke keep +whatever `--solver-name` selects — gurobi, cplex, xpress, glpk — and nothing here +constrains them. Each cylinder builds its own copy of the scenario models and solves +them with its own solver; the only coupling between cylinders is the numeric exchange +of `W` and bounds, which is solver-agnostic. + +The mechanism already exists. `apply_solver_specs(name, spoke, cfg)` overlays a +per-cylinder solver name onto the spoke's options dict: + +```python +# mpisppy/utils/cfg_vanilla.py +if _hasit(cfg, name+"_solver_name"): + options["solver_name"] = cfg.get(name+"_solver_name") +``` + +So the factory declares `ipopt_outer_bound_solver_name` with default `"ipopt"` and +calls `apply_solver_specs("ipopt_outer_bound", ...)` like every other spoke factory. +The guard in §6 rejects a value that does not name Ipopt. + +### 7.1 Global solver options must not leak into this cylinder + +This part is not free. `apply_solver_specs` deliberately layers per-spoke options *on +top of* the global `--solver-options` dict, and the global keys remain in place. That +is right for spokes sharing the global solver and wrong here, because **Ipopt +hard-fails on an unrecognized keyword** rather than ignoring it. Observed on this +branch: + +``` +ERROR: Solver log: Ipopt 3.13.2: ... Unknown keyword "acceptable_iter" +pyomo.common.errors.ApplicationError: Solver (ipopt) did not exit normally +``` + +So an entirely ordinary run — `--solver-name gurobi --solver-options "mipgap=0.01"` +for the hub, this cylinder attached alongside — would kill the cylinder on its first +solve, with an error naming Ipopt rather than the option routing. + +Decision: this cylinder does **not** inherit the global solver-options layer. It starts +from an empty base, and Ipopt-specific settings arrive through +`--ipopt-outer-bound-solver-options`. Filtering the global dict against a list of +Ipopt-known keywords was considered and rejected: the list would have to track Ipopt +releases, and silently dropping an option the user set is worse than never applying it. + +### 7.2 Convexity is a property of the model, not of the routing + +Worth stating because the flexibility above invites the wrong inference: the §6 model +guards apply to the shared scenario model, not to this cylinder's private copy of it. +If the model has integer variables this cylinder is inapplicable no matter what solver +anything else runs — and "another cylinder needs a MIP solver" is usually the signal +that it does. The routing lets Ipopt coexist with a MIP solver on a *convex* model (a +hub pushing an LP through gurobi, say); it does not let it certify a non-convex one. + +## 8. Combining across the cylinder's ranks + +This needs no new machinery, but it depends on an existing invariant that must not be +"optimized away", so it is recorded here. + +`SPOpt.Ebound()` already does exactly the required reduction — `Allreduce(..., MPI.SUM)` +of `p_s · outer_bound_s` over `self.mpicomm` — and already returns `None`, collectively, +if *any* scenario on *any* rank is missing its bound. + +**Every scenario in the sum must use the same `W`.** This is not fussiness. Let `W'` +mix generations across scenarios. At the true optimum `(x*, x̄*)`, + +``` + Σ_s p_s L_s(W'_s) ≤ Σ_s p_s [f_s(x*) + W'_sᵀ x*] = OPT + (Σ_s p_s W'_s)ᵀ x̄* +``` + +which bounds `OPT` only when `Σ_s p_s W'_s = 0` — a *joint* condition on the whole +weight vector that a mixture of generations does not satisfy. A cross-rank mixture +produces a number that is not a bound at all, and nothing downstream would detect it. + +Cross-rank agreement is already enforced: `get_receive_buffer(..., synchronize=True)` +routes through `_write_ids_agree`, which `Allreduce`s the write id and rejects a +mixed-generation read for retry. The new cylinder relies on this and must not pass +`synchronize=False`. + +The same argument forbids the stale-bound fallback used on the `outer_bound_only` path +(`spopt.py`: "Leave outer_bound at its previous value"), which is sound only if *every* +scenario is stale together. This cylinder instead sets `outer_bound = None` for any +scenario whose certificate fails, letting `Ebound()` return `None` and sending nothing. +Whether the existing Lagrangian spoke has the same exposure is a separate question, +noted in §12. + +## 9. Implementation sketch + +| File | Change | +|---|---| +| `mpisppy/utils/dual_certificate.py` | **Landed (Phase 1).** Pure function of a solved Pyomo model + its `dual` suffix → certified lower bound. No MPI, no cylinder, no PH. Named neutrally because only the §5 sign table is Ipopt-specific; the convention is a `sign_convention="ipopt"` argument rather than a hard-coded assumption. API: `check_model_is_certifiable`, `unbounded_variables`, `certified_lower_bound`, `CertificateError`. | +| `mpisppy/tests/test_dual_certificate.py` | **Landed (Phase 1).** Wired into `run_coverage.bash` and the `unit-tests` CI job. | +| `mpisppy/cylinders/ipopt_outer_bound.py` | **Landed.** `IpoptOuterBound(_LagrangianMixin, OuterBoundWSpoke)`; `outer_bound_only = False` (duals and primals are both needed); per-scenario certificate → `_mpisppy_data.outer_bound` → `Ebound()`. | +| `mpisppy/utils/config.py` | **Landed.** `ipopt_outer_bound_args()`. | +| `mpisppy/utils/cfg_vanilla.py` | **Landed.** `ipopt_outer_bound_spoke()` factory, alongside `lagrangian_spoke()`. | +| `mpisppy/generic/spokes.py` | **Landed.** Spoke registry. Note this is *not* `generic_cylinders.py`, where an earlier draft of this table put it — the driver delegates spoke construction to `generic/spokes.py`, and the arg registration to `generic/parsing.py`. | +| `mpisppy/generic/parsing.py` | **Landed.** Registers `ipopt_outer_bound_args()`. | +| `mpisppy/tests/test_ipopt_outer_bound.py` | **Landed.** Wiring tests (no solver, no MPI) plus an end-to-end run against the EF optimum. | +| `doc/src/spokes.rst` | **Landed.** User-facing page. | + +Prep attaches a `dual` Suffix (IMPORT) to each scenario. The objective expression after +`PH_Prep(attach_prox=False)` is already `f_s + W_sᵀx`, so `φ_s` builds directly on top +of it. + +The certificate engine is deliberately a standalone utility rather than a method on the +spoke: it is the part with the interesting math, and it is fully testable serially. +The *cylinder* keeps the Ipopt name, because the scope decision in §1 is real — only +Ipopt's conventions are measured and only Ipopt is accepted by the §6 guard. + +## 10. Phased rollout + +Each phase is its own review-sized PR and is green on its own. + +- **Phase 1 — certificate engine. DONE, in this design branch.** + `utils/dual_certificate.py` plus the §6 model guards plus + `tests/test_dual_certificate.py`. No cylinder, no MPI, no config surface. +- **Phase 2 — the cylinder. DONE.** `IpoptOuterBound`, config surface, `cfg_vanilla` + factory, driver wiring. +- **Phase 3 — parallel + docs. DONE.** Two-rank test exercising the `Ebound` + reduction; `doc/src/spokes.rst`; a driver command-line smoke run. + +**Phases 1–3 ship as a single PR**, revising the "each phase its own PR" plan above. +The reason is review latency rather than principle: with no reviewer available in this +area for some time, splitting buys nothing and costs coherence — and Phase 1 on its own +is arguably *harder* to review, being a mathematical argument with no consumer. Together, +a reviewer can watch the bound get produced against a known EF optimum. + +Deferred, and possibly forever — neither is needed for the feature to be complete, so +both are better filed as issues than kept as planned phases: + +- **Optional tightening re-solve.** Re-solve `min_{v∈B} φ_s(v)` when the correction + exceeds a threshold, and re-certify at the new point. By §3.4's own argument the + correction is ~0 at a converged solve, so this may never be worth building. +- **Maximization** through the concave mirror. A hard error today, which is a fine + permanent state unless someone actually needs it. + +### 10.1 Ipopt in CI — done, and it is the good build + +There was no Ipopt anywhere in `.github/workflows/`. Phase 1 was unaffected, since its +tests are deliberately solver-free (§11), but Phase 2 and Phase 3 could not have been +meaningfully green: their tests would have *skipped* rather than failed, which is the +worst of both. So this landed with Phase 1 rather than being left to Phase 3. + +**Job `ipopt-tests`**, in `test_pr_and_main.yml`. It pulls the release bundle from +`IDAES/idaes-ext` — the same source Pyomo's and Egret's own CI use, both of which were +read directly rather than reconstructed from memory: + +- Pyomo: `.github/workflows/test_pr_and_main.yml`, step "Install Ipopt" +- Egret: `.github/workflows/egret.yml`, step "Install Solvers" + +Why this source and not pip or conda: **pip/conda Ipopt is built against MUMPS only.** +The idaes-ext bundle additionally carries the HSL linear solvers, which is what makes +Ipopt usable on real NLPs. Measured on the extracted bundle: `ma27`, `ma57`, `ma97` all +solve; `ma86` does not ship. The job asserts those three actually work before running +any test, because a MUMPS-only fallback would otherwise show up much later as an +abnormal exit rather than as a clear "wrong build" message. + +Verified before committing, by rehearsing the whole step locally rather than trusting +the YAML: the tag extraction returns `3.4.2` from the live API and the hardcoded +fallback correctly triggers on a junk response; the tarball unpacks **flat** (so the +extraction directory itself goes on `PATH`); the extracted binary runs; `ma27/ma57/ma97` +all solve from it; and the Phase 1 suite passes against that exact binary. + +One non-obvious prerequisite: **the idaes-ext binary carries no RPATH** and resolves +`libgfortran`, `liblapack` and `libblas` from the system. A bare runner has none of +them, and the failure surfaces as an unhelpful load error at the first solve rather +than at install time, so the job `apt-get install`s the three before downloading. +Pyomo's workflow does the same, for the same reason. + +Two facts worth keeping: + +- `3.4.2` (2024-08-12) is still the current idaes-ext release, and is also what is + installed on the development machine — so local results and CI results are from the + same build, not merely from "some Ipopt". +- The bundle targets `ubuntu2204` while `ubuntu-latest` runners are 24.04. This is fine: + the same tarball is what runs on the 24.04 / glibc 2.39 development machine. + +## 11. Test plan + +**Phase 1's tests need no solver, by construction.** The certificate is a pure function +of the point the variables hold and the values in the `dual` suffix, so both are set by +hand and every expected number is exact analytic arithmetic rather than a recorded +observation. Given §10.1 this is not a stylistic preference: a solver-gated suite would +skip in CI and report nothing. The handful of assertions that genuinely require Ipopt — +that its reported dual *signs* are what the table assumes — are isolated in one +`skipUnless` class. + +The running example is `min (x−3)² + (y−2)² s.t. x + y ≤ 1, x,y ∈ [−10,10]`, with +optimum 8 at `(1, 0)` and multiplier 4, all analytic. + +- **Analytic**: `q̂ = 8` *exactly* at the KKT point — the correction term is 0 there, so + this is an equality assertion, not a tolerance. +- **Sign table**: one test per orientation (`<=`, `>=`, `==`, ranged), each with the + dual Ipopt reports for that orientation, all four required to return exactly 8. A lost + minus sign shows up as a wrong number, not a near-miss. +- **Degradation**: validity from points that are non-optimal, and from points that are + outright *infeasible* — the certificate requires neither. Under real truncated solves + (`max_iter ∈ {1,2,3,5,8}`, Ipopt-gated) validity is asserted at every level and the + converged bound is asserted to be at least as tight as each truncated one. Pairwise + monotonicity across `max_iter` is deliberately **not** asserted: the iterate path is a + solver detail that can differ by linear-solver build, and a test that fails on a + machine with HSL but passes on MUMPS is worse than no test. +- **Wrong-sign robustness**: feeding the `>=` dual to a `<=` constraint must yield a + *loose but valid* bound (−68 against an optimum of 8), pinning §5.1 down as a test + rather than a claim. +- **Box width**: widening the box from ±10 to ±1000 must loosen the bound and must not + invalidate it, pinning the §3.3 scaling law. +- **Guards**: integer var, binary var, nonlinear equality, maximize, no objective, + missing `dual` suffix, missing dual for a constraint, unknown sign convention — each + raises. Plus the negative cases that matter as much: a *fixed* discrete variable is + allowed (it is a constant, so it carries no convexity claim) and a *nonlinear + inequality* is allowed (convex inequalities are the entire point; only equalities must + be affine). The prox and non-Ipopt-solver guards are cylinder-level and land in + Phase 2. +- **Unbounded variable**: asserts the §6.1 path — `None` rather than `-inf`, and only + when the unbounded direction actually carries a nonzero gradient component. Includes + the half-open case, where whether the bound survives depends on which way the gradient + points. Also asserts the `fbbt` pre-pass rescues a variable whose bounds are implied + by the constraints. +- **Cushion**: `ε_rel = 0` reproduces `q̂` exactly; the default shaves it by + `1e-9·(1+|q̂|)` and never turns a valid bound invalid. +- **Non-finite results** (no solver): a NaN dual, a NaN in the point, and an infinite + dual must each yield `None` or a valid number, never a non-finite one. The `+inf` case + is the one that matters — NaN is rejected downstream by accident, `+inf` would be + latched as an improvement — so it is asserted explicitly rather than left to the + general check. +- **Ill-conditioning** (Ipopt-gated): Hilbert-matrix QPs at `cond(H)` ≈ 1.6e13 and + 3.5e17, with row scalings of 1e±8 and an objective offset of 1e12, certified at + `max_iter ∈ {1,2,3,5,10,∞}`. Each bound is checked against `f` at a point *constructed* + to be feasible by bisecting a uniform downshift — comparing against `f(v̂)` from the + solver instead is the trap described in §5.3 and produces false violations. Two + supporting assertions keep it honest: that the 1e−8 row scaling really does drive the + multiplier past 1e6 (or the cancellation variant would silently stop exercising + cancellation), and that the well-scaled case closes to 1e−6 relative while the badly + scaled one only has to stay valid — pinning §5.3's "tightness, not validity" as a test. + Tolerances are a few thousand ulps relative, which is the §5.2 arithmetic allowance; + an exact comparison at an offset of 1e12 would be testing the FPU. +- **Integration**: `farmer` is linear, hence convex, and Ipopt solves it; compare the + cylinder's bound against the EF optimum. +- **MPI**: 2-rank cylinder, same comparison, exercising `Ebound`'s reduction. The hub + runs Ipopt too by default so the test needs no MIP solver; a second case runs the hub + on a MIP solver and skips when none is available, which is what actually demonstrates + the §7 routing claim. +- **Option routing** (no solver, no MPI): that the global `--solver-options` layer does + *not* reach this spoke while the per-spoke layer does, and — the assertion that makes + that meaningful — that `lagrangian_spoke` still *does* inherit the global layer. This + is the §7.1 decision, and it is the part most likely to break silently. + +**No `run_all.py` entry.** The usual home for "the documented command line still works" +is `run_all.py`, but `do_one` has no skip machinery, so an entry there would force an +Ipopt install into both `run_all` CI jobs, which have no other use for one. The same +coverage is bought far more cheaply by a driver smoke run inside the `ipopt-tests` job, +which already has Ipopt. Worth revisiting if Ipopt ever lands in those jobs for another +reason. + +New `mpisppy/tests/test_*.py` files go into `run_coverage.bash` **and** +`.github/workflows/test_pr_and_main.yml` in the same commit, or codecov reports 0% on +the patch. + +## 12. Decisions + +Every question this design opened is now settled; item 4 turned out to be a bug in +shipped code rather than a design choice. + +1. **Cushion default — on.** `ε_rel = 1e-9`, `--ipopt-outer-bound-cushion` to change, + 0 to disable. See §5.2 for what it is and is not worth. +2. **Unbounded variables — warn, then report no bound.** `fbbt` first; anything still + unbounded produces a rank-0 setup warning naming the variables and a `None` bound, + not an exception. See §6.1. +3. **Naming — neutral engine, Ipopt cylinder.** `utils/dual_certificate.py` takes the + sign convention as an argument; `cylinders/ipopt_outer_bound.py` is the Ipopt-scoped + consumer. See §9. +4. **Pre-existing stale-bound exposure — confirmed, and fixed separately.** The + question was whether the existing stale-`outer_bound` fallback admits the §8 + mixed-`W` case in practice. It does. `solve_one` left a subproblem's previous + bound in place when a solve produced no bound, so `Ebound` could form a sum + mixing one scenario's stale bound with the others' fresh ones — not a bound at + all, and invisible to `Ebound`'s missing-bound check, which tests for `None` and + sees a number. The same exposure existed on the ordinary failed-solve path. + Reproduced and fixed in **PR #839**, based on `main` rather than stacked here, + since it is a soundness bug in shipped code and independent of this design. + Ipopt itself never reaches that path (it reports `Lower_bound = -inf`, not + `None`), so nothing in this design depends on the outcome. +5. **Non-finite results — reject, do not return.** A diverged solve can leave NaN in + the point or the duals, and an infinite multiplier yields NaN or `±inf`. All of it + used to propagate to the caller. The engine now tests `math.isfinite(q̂)` and returns + `None`, reusing the disposition it already has for an unbounded direction. The case + that makes this necessary rather than tidy is `+inf`: NaN is discarded downstream by + accident, since NaN loses every comparison and the hub's update test is `new > old`, + but `+inf` would compare as an *improvement* and be latched as an outer bound that is + not one. See §5.2. +6. **Box tightening — two mechanisms, two arguments.** `fbbt` at setup and the + nonant-bounds channel each iteration both shrink the box the certificate minimizes + over, and they are not sound for the same reason. §3.6 states the condition each has + to meet and which one meets which. Recorded as a decision because the tempting + one-line justification is wrong in a way that would not show up in any test. diff --git a/doc/src/spokes.rst b/doc/src/spokes.rst index 532bde125..d70b25c09 100644 --- a/doc/src/spokes.rst +++ b/doc/src/spokes.rst @@ -34,6 +34,120 @@ hedging algorithm for stochastic mixed-integer programs` by Gade et al [gade2016]_. It takes W values from the hub and uses them to compute a bound. +ipopt_outer_bound +^^^^^^^^^^^^^^^^^ + +An outer bound for problems whose scenario subproblems are **convex NLPs solved +with Ipopt**, enabled with ``--ipopt-outer-bound``. + +The Lagrangian spoke gets its bound from the solver's dual bound. Ipopt is not a +branch-and-bound solver and reports none, so on a convex NLP that spoke produces +nothing usable. This one computes the bound itself, from the subproblem's own +duals. + +This spoke does not simply report the solved objective value. That value is +measured *at a point*, hence an inner bound for a minimization -- the wrong +direction. The spoke instead computes a Lagrangian weak-duality bound, corrected +by a tangent-plane underestimator minimized in closed form over the variable box. +The result is valid for *any* multipliers, so no assumption that Ipopt converged +is needed: a truncated or sloppy solve gives a loose bound rather than a wrong +bound. At an exact KKT point the correction vanishes and the bound equals the +subproblem optimum. + +Cost is one solve per scenario per iteration, the same as the Lagrangian spoke. + +.. warning:: + **Convexity is assumed and mostly cannot be checked.** If the objective is + non-convex, or an inequality is non-convex *in canonical form*, the bound is + simply wrong rather than merely loose. + + Canonical form matters, and it is easy to get backwards. Every inequality is + rewritten as ``g(v) <= 0``, which negates the body of a ``>=`` row: + + ========================== ==================== ========================== + as written canonical ``g`` requirement on the body + ========================== ==================== ========================== + ``body <= upper`` ``body - upper`` convex + ``body >= lower`` ``lower - body`` **concave** + ``lo <= body <= up`` both of the above affine + ``body == rhs`` ``body - rhs`` affine + ========================== ==================== ========================== + + So the theorem applies to ``x**2 <= 4`` but not to ``x**2 >= 1``, even though + both are written with a convex body -- and the feasible set of the latter is + not convex at all. Getting this wrong yields an outer bound that can exceed + the true optimum. + + What *is* checked, as a hard error at setup: discrete variables, nonlinear + equality constraints, nonlinear two-sided (ranged) constraints, a + maximization objective, a solver that is not Ipopt, and a ``dual`` Suffix + the scenario creator already attached in a direction that does not import + (the certificate needs the solver's duals back, so ``Suffix.IMPORT`` or + ``Suffix.IMPORT_EXPORT`` is required). The affine cases are decidable, so + they are enforced; convexity of a one-sided nonlinear body is not, so + convexity remains a user assertion. + +Two things determine whether the bound is any good: + +**Variable bounds.** The certificate minimizes over the box of variable bounds, +so its looseness is roughly ``sum_i |d_i phi| * (width of the box in the +descending direction)``. Tight bounds give a tight bound; enormous ones give a +valid but useless number. ``fbbt`` runs first to recover bounds implied by the +constraints. A variable still unbounded afterwards produces a warning at setup, +and the spoke reports nothing on iterations where that variable's gradient +component is nonzero -- it stays quiet rather than sending a wrong number. + +**Solver options.** Unlike every other spoke, this one does **not** inherit the +global ``--solver-options``. Ipopt hard-fails on an unrecognized keyword rather +than ignoring it, so a perfectly ordinary run (a MIP solver and its options for +the hub, this spoke attached alongside) would otherwise kill the spoke on its +first solve. Pass Ipopt settings through ``--ipopt-outer-bound-solver-options``. + +The hub and the other spokes keep whatever ``--solver-name`` selects; only this +spoke is pinned to Ipopt, via ``--ipopt-outer-bound-solver-name`` (default +``ipopt``). Note that this routing lets Ipopt coexist with a MIP solver on a +*convex* model -- it does not let it certify a non-convex one. If the model has +integer variables this spoke is inapplicable no matter what anything else runs. + +**An inexact or ill-conditioned solve is safe.** The certificate assumes nothing +about the accuracy of the solve. The model is convex by assumption, so it has no +non-global local minima, and Ipopt returning a sub-optimal answer can only mean +it stopped short of converging -- an inexact point with inexact multipliers, both +of which the underlying theorem admits. Ill-conditioning does not enter either: +the certificate evaluates the objective and one gradient, and inverts nothing, so +there is no linear solve for a condition number to amplify. What both cost is +tightness. The looseness term grows as the point moves away from optimal, so a +badly conditioned subproblem reports a weak bound, and the hub keeps the best +outer bound it has seen and ignores it. + +``--ipopt-outer-bound-cushion`` (default ``1e-9``) subtracts a small relative +amount, ``q - eps*(1+|q|)``, from the reported bound. This is last-bit hygiene +against floating point, not a proof-carrying margin; pass ``0`` to disable it. + +The one case worth raising it for is a badly *scaled* model. The bound is summed +as ``f + lam^T g + mu^T h``, so its rounding error tracks the size of those terms, +while the cushion tracks the size of the answer. Multipliers of order ``1e8`` +against an objective of order one put the error floor near ``1e-8`` while the +default cushion is ``1e-9``. Rescaling the offending constraint rows is the better +fix; raising the cushion is the cheap one. A solve that diverges outright produces +no number at all -- a non-finite result is rejected and the spoke stays quiet. + +Maximization is not supported and raises at setup. + +.. note:: + Ipopt builds obtained from the IDAES ``idaes-ext`` distribution (the usual way + to get one with good linear solvers, and what mpi-sppy's CI installs) link the + Harwell Subroutine Library and **default to the** ``ma27`` **linear solver** + rather than to MUMPS. That is worth knowing because results can differ + slightly between them, and because HSL asks that its use be acknowledged: + + HSL, a collection of Fortran codes for large-scale scientific computation. + See https://www.hsl.rl.ac.uk/ + + Pass ``--ipopt-outer-bound-solver-options "linear_solver=mumps"`` to choose + otherwise. + + Subgradient ^^^^^^^^^^^ diff --git a/mpisppy/cylinders/ipopt_outer_bound.py b/mpisppy/cylinders/ipopt_outer_bound.py new file mode 100644 index 000000000..9cd0961b5 --- /dev/null +++ b/mpisppy/cylinders/ipopt_outer_bound.py @@ -0,0 +1,276 @@ +############################################################################### +# mpi-sppy: MPI-based Stochastic Programming in PYthon +# +# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for +# Sustainable Energy, LLC, The Regents of the University of California, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for +# full copyright and license information. +############################################################################### +# An outer-bound spoke for convex NLP subproblems solved with Ipopt. +# +# The ordinary Lagrangian spoke reads its bound off the solver's dual bound. +# Ipopt is not a branch-and-bound solver and reports none, so that spoke warns +# and produces nothing usable. This one computes the bound itself, from the +# subproblem's own duals, using mpisppy.utils.dual_certificate -- see that +# module for the argument. The short version: the value at the returned point +# is an *inner* bound, and what makes an outer bound available without any +# convergence assumption is Lagrangian weak duality plus a tangent-plane +# underestimator minimized in closed form over the variable box. +# +# Convexity of the scenario subproblems is the user's assertion. The parts of it +# that can be checked mechanically are checked at setup and are hard errors; see +# _check_setup_guards. + +import warnings + +import pyomo.environ as pyo + +import mpisppy.utils.sputils as sputils +from mpisppy.cylinders.lagrangian_bounder import _LagrangianMixin +from mpisppy.cylinders.spoke import OuterBoundWSpoke +from mpisppy.utils.dual_certificate import ( + CertificateError, + certified_lower_bound, + check_model_is_certifiable, + unbounded_variables, +) + + +class IpoptOuterBound(_LagrangianMixin, OuterBoundWSpoke): + + # 'N' for NLP. Not 'I': that is InnerBoundSpoke's character, and the hub + # prints the outer and inner chars side by side, so an 'I' in the outer + # column would read as an inner bound. + converger_spoke_char = 'N' + + # The certificate reads the point *and* the duals off the solved model, so + # unlike the Lagrangian spoke this one cannot skip loading the solution. + outer_bound_only = False + + def ipopt_outer_bound_prep(self): + """lagrangian_prep plus what the certificate needs: a dual suffix on + every subproblem, the setup guards, and the bound tightening.""" + self.lagrangian_prep() + + # PH_Prep(attach_prox=False) is what makes this the Lagrangian + # relaxation rather than a proximal subproblem, so check that the prox + # term is not switched on. prox_on is created by attach_Ws_and_prox and + # starts at 0, so this normally passes; it is here to catch a future + # path that enables it, in which case the number below would not be a + # Lagrangian bound at all. + first = next(iter(self.opt.local_scenarios.values()), None) + if (first is not None + and hasattr(first._mpisppy_model, "prox_on") + and not self.opt.prox_disabled): + raise CertificateError( + "ipopt_outer_bound requires the proximal term to be off; " + "the bound it computes is a Lagrangian bound and a proximal " + "subproblem is not the Lagrangian relaxation" + ) + + for s in self.opt.local_scenarios.values(): + # Existence is not enough: a scenario_creator may already attach an + # EXPORT or LOCAL `dual` suffix (a common way to supply dual warm + # starts), and reusing that would import nothing. + existing = getattr(s, "dual", None) + if existing is None: + s.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT) + elif not existing.import_enabled(): + raise CertificateError( + f"scenario {s.name} already has a `dual` Suffix that does " + "not import; the certificate needs the solver's duals. Use " + "Suffix.IMPORT or Suffix.IMPORT_EXPORT." + ) + + self._warned = set() + self._check_setup_guards() + + # Snapshot which nonants are fixed now, so a fixing extension that + # fixes more of them later can be caught: fixing a nonant restricts the + # subproblem, which can only raise its minimum, so the result would + # bound the restricted problem and not the original. Nonants the + # scenario creator fixed are part of the problem and are fine. + self._fixed_at_setup = { + (sname, ndn_i): xvar.fixed + for sname, s in self.opt.local_scenarios.items() + for ndn_i, xvar in s._mpisppy_data.nonant_indices.items() + } + + def _check_setup_guards(self): + """Hard errors for the parts of the theorem that are checkable, and a + warning for the part that only costs tightness.""" + solver_name = self.opt.options.get("solver_name") or "" + if "ipopt" not in solver_name: + raise CertificateError( + f"ipopt_outer_bound is scoped to Ipopt, but its solver is " + f"{solver_name!r}. The dual sign conventions it relies on are " + "measured from Ipopt only. Set --ipopt-outer-bound-solver-name." + ) + + for sname, s in self.opt.local_scenarios.items(): + try: + check_model_is_certifiable(s) + except CertificateError as e: + raise CertificateError(f"scenario {sname}: {e}") from None + + # fbbt first (it can only shrink the box, which makes the bound + # tighter without dropping a feasible point), then say something if a + # variable is still unbounded. Not an error: this cylinder is an + # optional source of a bound, and a model that is merely under-bounded + # is not a broken model. It may simply report nothing. + still_unbounded = {} + for sname, s in self.opt.local_scenarios.items(): + names = unbounded_variables(s, do_fbbt=True) + if names: + still_unbounded[sname] = names + if still_unbounded and self.cylinder_rank == 0: + sname, names = next(iter(still_unbounded.items())) + warnings.warn( + f"ipopt_outer_bound: {len(still_unbounded)} scenario(s) have " + "variables with no finite bound after fbbt, for example " + f"{sname}: {', '.join(names[:5])}" + f"{' ...' if len(names) > 5 else ''}. The certificate minimizes " + "over the variable box, so an unbounded direction with a " + "nonzero gradient yields no bound and this spoke will stay " + "quiet on those iterations. Bounding those variables is what " + "makes this spoke useful." + ) + + def _warn_once(self, key, message): + """Rank-0, once-per-run warning. This spoke is an optional source of a + bound, so it complains and stands down rather than taking the run with + it -- an exception here would propagate out of the iteration loop and + MPI_Abort the hub and every other spoke.""" + if key in self._warned: + return + self._warned.add(key) + if self.cylinder_rank == 0: + warnings.warn(message) + + def _nonants_newly_fixed(self): + """True if anything fixed a nonant since setup, in which case no bound + can be reported: fixing restricts the subproblem, so its minimum bounds + the restricted problem and not the original.""" + newly = [ + f"{sname}:{ndn_i}" + for sname, s in self.opt.local_scenarios.items() + for ndn_i, xvar in s._mpisppy_data.nonant_indices.items() + if xvar.fixed and not self._fixed_at_setup[(sname, ndn_i)] + ] + if not newly: + return False + self._warn_once( + "fixed_nonants", + "ipopt_outer_bound: nonanticipative variables were fixed after " + f"setup ({', '.join(newly[:5])}" + f"{' ...' if len(newly) > 5 else ''}). Fixing restricts the " + "subproblem, so its minimum is a bound on the restricted problem " + "and not on the original. This spoke will report no bound until " + "they are unfixed; remove the fixing extension from this spoke." + ) + return True + + def _solve_and_certify(self, warmstart=sputils.WarmstartStatus.PRIOR_SOLUTION): + """Solve every subproblem, then replace the solver's (useless) bound + with the certificate. Returns the expected outer bound, or None.""" + # This shrinks the box the certificate minimizes over, so it needs an + # argument. Note the tempting one -- "a smaller box removes points, and + # fewer points can only raise an infimum" -- is an argument that the + # bound gets TIGHTER, which is precisely the direction that could break + # it. What is actually needed is that every scenario's box still holds + # one COMMON optimal solution x* of the full problem: then the + # certificate is below phi_s(x*) <= f_s(x*) + W_s'x* for each s, and the + # p-weighted sum is below OPT. The fbbt done at setup gives this the + # easy way, by removing no feasible point at all. This channel does not + # -- only reduced_costs_spoke sends it, and reduced-cost fixing does + # discard feasible points -- so it rests on that spoke's own contract + # that an optimal solution survives, plus the fact that the bounds are + # broadcast and applied identically to every scenario, which is what + # makes the surviving solution common rather than per-scenario. A + # sender that guarantees neither would break the bound silently. + self.receive_nonant_bounds() + verbose = self.opt.options['verbose'] + teeme = self.opt.options.get('tee-rank0-solves', False) + + self.opt.solve_loop( + solver_options=self.opt._effective_solver_options(self.opt._PHIter), + dtiming=False, + gripe=True, + tee=teeme, + verbose=verbose, + need_solution=True, + warmstart=warmstart, + ) + + if self._nonants_newly_fixed(): + for s in self.opt.local_scenarios.values(): + s._mpisppy_data.outer_bound = None + return self.opt.Ebound(verbose) + + for s in self.opt.local_scenarios.values(): + # solve_loop has just written results.Problem[0].Lower_bound here, + # which for Ipopt is -inf. Overwrite it with the certificate, or + # with None when there is no certificate to be had -- Ebound then + # declines collectively rather than folding a -inf into the sum. + if not s._mpisppy_data.solution_available: + s._mpisppy_data.outer_bound = None + continue + try: + s._mpisppy_data.outer_bound = certified_lower_bound( + s, sign_convention="ipopt", eps_rel=self._cushion) + except CertificateError as e: + # Routine solver outcomes can leave a constraint without a dual. + # Report no bound, as the module's contract says, rather than + # taking down the hub and every other spoke from inside the + # iteration loop. + self._warn_once( + "certificate_failed", + f"ipopt_outer_bound: no certificate for {s.name} ({e}); " + "reporting no bound this iteration.") + s._mpisppy_data.outer_bound = None + + return self.opt.Ebound(verbose) + + @property + def _cushion(self): + return self.opt.options.get("ipopt_outer_bound_cushion", 1e-9) + + def _set_weights_and_solve(self, warmstart=sputils.WarmstartStatus.PRIOR_SOLUTION): + self.opt.W_from_flat_list(self.localWs) + return self._solve_and_certify(warmstart=warmstart) + + def main(self): + self.verbose = self.opt.options['verbose'] + extensions = self.opt.extensions is not None + + self.ipopt_outer_bound_prep() + + if extensions: + self.opt.extobject.pre_iter0() + + self.opt._PHIter = 0 + self.trivial_bound = self._solve_and_certify( + warmstart=sputils.WarmstartStatus.USER_SOLUTION) + + if extensions: + self.opt.extobject.post_iter0() + self.opt._PHIter += 1 + self.opt.current_solver_options = {} + + if self.trivial_bound is not None: + self.send_bound(self.trivial_bound) + if extensions: + self.opt.extobject.post_iter0_after_sync() + + while not self.got_kill_signal(): + if self.update_Ws(): + if extensions: + self.opt.extobject.miditer() + bound = self._set_weights_and_solve() + if extensions: + self.opt.extobject.enditer() + if bound is not None: + self.send_bound(bound) + if extensions: + self.opt.extobject.enditer_after_sync() + self.opt._PHIter += 1 diff --git a/mpisppy/generic/admm.py b/mpisppy/generic/admm.py index b62d065d8..226f9b4a0 100644 --- a/mpisppy/generic/admm.py +++ b/mpisppy/generic/admm.py @@ -56,7 +56,7 @@ def _count_cylinders(cfg): """ count = 1 # the hub spoke_flags = [ - "fwph", "lagrangian", "ph_dual", "relaxed_ph", + "fwph", "lagrangian", "ipopt_outer_bound", "ph_dual", "relaxed_ph", "subgradient", "xhatshuffle", "xhatxbar", "reduced_costs", ] for flag in spoke_flags: diff --git a/mpisppy/generic/parsing.py b/mpisppy/generic/parsing.py index a46c7cf0b..d76ebd099 100644 --- a/mpisppy/generic/parsing.py +++ b/mpisppy/generic/parsing.py @@ -126,6 +126,7 @@ def add_decomp_args(cfg): cfg.ph_xfeas_spoke_args() cfg.fwph_args() cfg.lagrangian_args() + cfg.ipopt_outer_bound_args() cfg.subgradient_bounder_args() cfg.xhatshuffle_args() cfg.xhatxbar_args() diff --git a/mpisppy/generic/spokes.py b/mpisppy/generic/spokes.py index d86da280a..a32d1ada5 100644 --- a/mpisppy/generic/spokes.py +++ b/mpisppy/generic/spokes.py @@ -58,6 +58,15 @@ def build_spoke_list(cfg, beans, scenario_creator_kwargs, vanilla.add_gapper(lagrangian_spoke, cfg, "lagrangian") lagrangian_spoke["rank_ratio"] = cfg.lagrangian_rank_ratio + # Certified outer bound for convex NLP subproblems solved with Ipopt + if cfg.ipopt_outer_bound: + ipopt_outer_bound_spoke = vanilla.ipopt_outer_bound_spoke(*beans, + scenario_creator_kwargs=scenario_creator_kwargs, + rho_setter=rho_setter, + all_nodenames=all_nodenames, + ) + ipopt_outer_bound_spoke["rank_ratio"] = cfg.ipopt_outer_bound_rank_ratio + # dual ph spoke if cfg.ph_dual: ph_dual_spoke = vanilla.ph_dual_spoke(*beans, @@ -173,6 +182,8 @@ def build_spoke_list(cfg, beans, scenario_creator_kwargs, list_of_spoke_dict.append(fw_spoke) if cfg.lagrangian: list_of_spoke_dict.append(lagrangian_spoke) + if cfg.ipopt_outer_bound: + list_of_spoke_dict.append(ipopt_outer_bound_spoke) if cfg.ph_dual: list_of_spoke_dict.append(ph_dual_spoke) if cfg.relaxed_ph: diff --git a/mpisppy/tests/test_dual_certificate.py b/mpisppy/tests/test_dual_certificate.py new file mode 100644 index 000000000..8b83f70db --- /dev/null +++ b/mpisppy/tests/test_dual_certificate.py @@ -0,0 +1,556 @@ +############################################################################### +# mpi-sppy: MPI-based Stochastic Programming in PYthon +# +# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for +# Sustainable Energy, LLC, The Regents of the University of California, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for +# full copyright and license information. +############################################################################### +# Tests for mpisppy/utils/dual_certificate.py +# +# Most of these need no solver. The certificate is a pure function of the point +# the variables hold and the values in the `dual` suffix, so setting both by +# hand gives exact arithmetic and covers the module without depending on ipopt +# being installed -- which it is not, in CI. The solver-dependent tests are +# gathered in TestWithIpopt and skip cleanly. +# +# The running example throughout is +# +# min (x-3)^2 + (y-2)^2 s.t. x + y <= 1, x, y in [-10, 10] +# +# whose optimum is x=1, y=0 with value 8 and multiplier 4 -- all analytic, so +# every expected number below is exact rather than a recorded observation. + +import math +import unittest + +import pyomo.environ as pyo + +from mpisppy.tests.utils import announce_hsl_if_used +from mpisppy.utils.dual_certificate import ( + CertificateError, + certified_lower_bound, + check_model_is_certifiable, + unbounded_variables, +) + +OPT = 8.0 # analytic optimum of the running example +MULTIPLIER = 4.0 # analytic multiplier of the active constraint + +ipopt_available = pyo.SolverFactory("ipopt").available(exception_flag=False) + +if ipopt_available: + announce_hsl_if_used() + + +def _model(kind="le", bounds=(-10, 10), y_bounds=None): + """The running example, with the single constraint written four ways. + + All four describe the same feasible set, so all four have the same optimum + and the same |multiplier| -- only the sign Pyomo reports differs. + """ + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=bounds, initialize=0.0) + m.y = pyo.Var(bounds=y_bounds if y_bounds is not None else bounds, initialize=0.0) + m.obj = pyo.Objective(expr=(m.x - 3) ** 2 + (m.y - 2) ** 2, sense=pyo.minimize) + if kind == "le": + m.c = pyo.Constraint(expr=m.x + m.y <= 1) + elif kind == "ge": + m.c = pyo.Constraint(expr=-m.x - m.y >= -1) + elif kind == "eq": + m.c = pyo.Constraint(expr=m.x + m.y == 1) + elif kind == "range": + m.c = pyo.Constraint(expr=pyo.inequality(-10, m.x + m.y, 1)) + else: + raise ValueError(kind) + m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT) + return m + + +def _place(m, x, y, dual): + """Put the model at a point with a given dual, as a solve would have.""" + m.x.value = x + m.y.value = y + m.dual[m.c] = dual + return m + + +# The dual value ipopt reports for each orientation when the true multiplier is +# +4. Verified against ipopt in TestWithIpopt.test_sign_table_matches_ipopt. +IPOPT_DUAL_AT_OPTIMUM = {"le": -4.0, "ge": 4.0, "eq": -4.0, "range": -4.0} + + +class TestCertificateMath(unittest.TestCase): + """No solver required: the point and the duals are supplied directly.""" + + def test_exact_at_kkt_point(self): + # At the exact KKT point the correction term is identically zero and the + # certificate collapses to the optimal value -- strong duality, exactly. + m = _place(_model("le"), 1.0, 0.0, IPOPT_DUAL_AT_OPTIMUM["le"]) + self.assertEqual(certified_lower_bound(m, eps_rel=0.0), OPT) + + def test_sign_table_all_orientations(self): + # Each orientation needs a different rule to recover the same canonical + # multiplier. A lost minus sign here shows up as a loose bound, so the + # assertion is equality, not just validity. + for kind, dual in IPOPT_DUAL_AT_OPTIMUM.items(): + with self.subTest(kind=kind): + m = _place(_model(kind), 1.0, 0.0, dual) + self.assertEqual(certified_lower_bound(m, eps_rel=0.0), OPT) + + def test_inactive_constraint_drops_out(self): + # Complementarity: a zero dual must contribute nothing. Unconstrained + # optimum is x=3, y=2 with value 0. + m = _model("le") + m.c.deactivate() + m.d = pyo.Constraint(expr=m.x + m.y <= 100) + m.dual[m.d] = 0.0 + _place_x, _place_y = 3.0, 2.0 + m.x.value, m.y.value = _place_x, _place_y + self.assertEqual(certified_lower_bound(m, eps_rel=0.0), 0.0) + + def test_valid_at_a_point_that_is_not_optimal(self): + # phi(0,0) = 9 + 4 + 4*(0+0-1) = 9; grad = (-2, 0); the box term drives + # x to its upper bound 10, contributing -2*(10-0) = -20. + m = _place(_model("le"), 0.0, 0.0, IPOPT_DUAL_AT_OPTIMUM["le"]) + q = certified_lower_bound(m, eps_rel=0.0) + self.assertEqual(q, -11.0) + self.assertLessEqual(q, OPT) + + def test_infeasible_point_still_gives_a_valid_bound(self): + # vhat need not be feasible: x+y = 4 violates x+y <= 1. The certificate + # does not care, which is why a truncated solve is safe. + m = _place(_model("le"), 2.0, 2.0, IPOPT_DUAL_AT_OPTIMUM["le"]) + self.assertLessEqual(certified_lower_bound(m, eps_rel=0.0), OPT) + + def test_wrong_sign_is_loose_never_wrong(self): + # Feed the `>=` dual to a `<=` constraint. max(-(+4), 0) = 0 clips the + # multiplier away, so phi degenerates to f -- loose, still valid. This + # is the robustness property that makes the whole approach safe. + m = _place(_model("le"), 1.0, 0.0, -IPOPT_DUAL_AT_OPTIMUM["le"]) + q = certified_lower_bound(m, eps_rel=0.0) + self.assertLess(q, OPT) + self.assertEqual(q, -68.0) + + def test_looseness_scales_with_box_width(self): + # The correction is |grad| times the distance to the far end of the box, + # so wide bounds give a valid but weak bound. Practical consequence: + # tight variable bounds are what make this cylinder useful. + narrow = certified_lower_bound( + _place(_model("le", bounds=(-10, 10)), 0.0, 0.0, -4.0), eps_rel=0.0 + ) + wide = certified_lower_bound( + _place(_model("le", bounds=(-1000, 1000)), 0.0, 0.0, -4.0), eps_rel=0.0 + ) + self.assertLess(wide, narrow) + self.assertLessEqual(narrow, OPT) + self.assertLessEqual(wide, OPT) + + def test_unbounded_direction_returns_none_not_minus_inf(self): + # d(phi)/dy = 2*(y-2) + lam = 2y when lam = 4, so any y != 0 gives y a + # nonzero gradient component. With y unbounded the box minimization is + # then -inf, and the contract is None ("no bound this time") rather than + # a -inf that would poison Ebound's sum. + m = _place(_model("le", y_bounds=(None, None)), 0.0, 1.0, -4.0) + self.assertIsNone(certified_lower_bound(m)) + + def test_unbounded_on_the_unused_side_only_still_bounds(self): + # Half-open is enough when the gradient points the other way: here + # d(phi)/dy = 2 > 0, so only the lower bound is consulted. + m = _place(_model("le", y_bounds=(-10, None)), 0.0, 1.0, -4.0) + self.assertIsNotNone(certified_lower_bound(m)) + m2 = _place(_model("le", y_bounds=(None, 10)), 0.0, 1.0, -4.0) + self.assertIsNone(certified_lower_bound(m2)) + + def test_unbounded_direction_with_zero_gradient_still_bounds(self): + # An infinite bound only defeats the certificate if that variable's + # gradient component is nonzero. At the KKT point it is exactly zero + # here, so a bound is still available. + m = _place(_model("le", y_bounds=(None, None)), 1.0, 0.0, -4.0) + self.assertEqual(certified_lower_bound(m, eps_rel=0.0), OPT) + + def test_cushion(self): + m = _place(_model("le"), 1.0, 0.0, -4.0) + exact = certified_lower_bound(m, eps_rel=0.0) + cushioned = certified_lower_bound(m) # default 1e-9 + self.assertEqual(exact, OPT) + self.assertLess(cushioned, exact) + self.assertAlmostEqual(exact - cushioned, 1e-9 * (1.0 + OPT), places=15) + + def test_nan_dual_yields_no_bound(self): + # A diverged solve can leave NaN in the `dual` suffix, and NaN + # propagates silently through phi and the correction. The contract is + # None, the same word this module already uses for "no bound this time". + m = _place(_model("le"), 1.0, 0.0, float("nan")) + self.assertIsNone(certified_lower_bound(m, eps_rel=0.0)) + + def test_nan_in_the_point_yields_no_bound(self): + m = _place(_model("le"), 1.0, 0.0, -4.0) + m.x.set_value(float("nan"), skip_validation=True) + self.assertIsNone(certified_lower_bound(m, eps_rel=0.0)) + + def test_infinite_duals_yield_no_bound_or_a_valid_one(self): + # -inf clips to lam = +inf, and inf*0 is NaN; +inf clips to lam = 0, + # which just drops the constraint and leaves a loose but valid bound. + # Neither may return a number above the optimum. + for d in (float("inf"), float("-inf")): + with self.subTest(dual=d): + q = certified_lower_bound(_place(_model("le"), 1.0, 0.0, d), + eps_rel=0.0) + self.assertTrue(q is None or q <= OPT) + + def test_a_non_finite_bound_never_reaches_the_caller(self): + # The guard is on the result, so it catches every route to a non-finite + # qhat, not just the ones enumerated above. +inf matters most: unlike + # NaN it would compare as an *improvement* to a hub tracking the best + # outer bound seen. + for x, y, d in ((float("nan"), 0.0, -4.0), + (1.0, float("nan"), -4.0), + (1.0, 0.0, float("-inf")), + (float("inf"), 0.0, -4.0)): + with self.subTest(x=x, y=y, dual=d): + m = _model("le") + m.x.set_value(x, skip_validation=True) + m.y.set_value(y, skip_validation=True) + m.dual[m.c] = d + q = certified_lower_bound(m, eps_rel=0.0) + self.assertTrue(q is None or math.isfinite(q)) + + def test_fixed_variables_are_constants(self): + # A fixed variable is not free in the box and must not contribute a + # correction term. Fix y at its optimal value; the bound is unchanged. + m = _place(_model("le"), 1.0, 0.0, -4.0) + m.y.fix(0.0) + self.assertEqual(certified_lower_bound(m, eps_rel=0.0), OPT) + + +class TestUnboundedVariables(unittest.TestCase): + def test_reports_variables_without_finite_bounds(self): + m = _model("le", y_bounds=(None, None)) + self.assertEqual(unbounded_variables(m, do_fbbt=False), ["y"]) + + def test_fbbt_rescues_bounds_implied_by_constraints(self): + m = _model("le", y_bounds=(None, None)) + m.b1 = pyo.Constraint(expr=m.y >= 0) + m.b2 = pyo.Constraint(expr=m.y <= 5) + self.assertEqual(unbounded_variables(m, do_fbbt=False), ["y"]) + self.assertEqual(unbounded_variables(m, do_fbbt=True), []) + self.assertEqual((m.y.lb, m.y.ub), (0, 5)) + + def test_fully_bounded_model_reports_nothing(self): + self.assertEqual(unbounded_variables(_model("le")), []) + + +class TestGuards(unittest.TestCase): + def test_clean_model_passes(self): + check_model_is_certifiable(_model("le")) # must not raise + + def test_integer_variable(self): + m = _model("le") + m.z = pyo.Var(bounds=(0, 10), domain=pyo.Integers) + with self.assertRaisesRegex(CertificateError, "discrete"): + check_model_is_certifiable(m) + + def test_binary_variable(self): + m = _model("le") + m.z = pyo.Var(domain=pyo.Binary) + with self.assertRaisesRegex(CertificateError, "discrete"): + check_model_is_certifiable(m) + + def test_fixed_discrete_variable_is_allowed(self): + # Fixed means constant, so it carries no convexity claim. + m = _model("le") + m.z = pyo.Var(bounds=(0, 10), domain=pyo.Integers) + m.z.fix(3) + check_model_is_certifiable(m) # must not raise + + def test_nonlinear_equality(self): + m = _model("le") + m.bad = pyo.Constraint(expr=m.x * m.y == 1) + with self.assertRaisesRegex(CertificateError, "nonlinear equality"): + check_model_is_certifiable(m) + + def test_nonlinear_le_inequality_is_allowed(self): + # Convex inequalities are the point of the exercise. For a `<=` row the + # canonical g is body - upper, so a convex body is what the theorem + # wants -- and convexity is the caller's assertion, not checkable here. + m = _model("le") + m.ok = pyo.Constraint(expr=m.x**2 + m.y**2 <= 100) + check_model_is_certifiable(m) # must not raise + + def test_nonlinear_ge_inequality_is_allowed_but_needs_a_concave_body(self): + # The sign trap. For a `>=` row the canonical g is lower - body, so the + # BODY must be concave, not convex. That is still the caller's + # assertion, so the guard lets it through -- but a convex body here + # (x**2 >= 1, which looks perfectly ordinary) makes phi non-convex and + # the certificate can then exceed the true optimum. The rule is stated + # in the module docstring and in spokes.rst; this test exists to pin + # down that the guard does NOT claim to catch it. + m = _model("le") + m.ok = pyo.Constraint(expr=m.x**2 >= 1) + check_model_is_certifiable(m) # must not raise + + def test_nonlinear_ranged_constraint_is_rejected(self): + # A two-sided row splits into both g = body - upper and + # g = lower - body, so the body would have to be convex AND concave -- + # affine. Unlike one-sided convexity, that IS decidable, so it is a + # hard error rather than an assertion. + m = _model("le") + m.bad = pyo.Constraint(expr=pyo.inequality(1, m.x**2, 100)) + with self.assertRaisesRegex(CertificateError, "two-sided"): + check_model_is_certifiable(m) + + def test_affine_ranged_constraint_is_allowed(self): + m = _model("le") + m.ok = pyo.Constraint(expr=pyo.inequality(1, 2 * m.x + m.y, 100)) + check_model_is_certifiable(m) # must not raise + + def test_maximize(self): + m = _model("le") + m.obj.sense = pyo.maximize + with self.assertRaisesRegex(CertificateError, "minimize-only"): + check_model_is_certifiable(m) + + def test_no_objective(self): + m = _model("le") + m.obj.deactivate() + with self.assertRaisesRegex(CertificateError, "exactly one active Objective"): + check_model_is_certifiable(m) + + def test_missing_dual_suffix(self): + m = _model("le") + m.del_component(m.dual) + with self.assertRaisesRegex(CertificateError, "no `dual` Suffix"): + certified_lower_bound(m) + + def test_missing_dual_for_a_constraint(self): + m = _model("le") # suffix present but never populated + m.x.value, m.y.value = 1.0, 0.0 + with self.assertRaisesRegex(CertificateError, "no dual available"): + certified_lower_bound(m) + + def test_unknown_sign_convention(self): + m = _place(_model("le"), 1.0, 0.0, -4.0) + with self.assertRaisesRegex(CertificateError, "unknown sign convention"): + certified_lower_bound(m, sign_convention="no-such-solver") + + +@unittest.skipUnless(ipopt_available, "ipopt is not available") +class TestWithIpopt(unittest.TestCase): + """The parts that can only be checked against the real solver: that ipopt's + reported dual signs are what the table in dual_certificate.py assumes.""" + + @staticmethod + def _solve(m, max_iter=None): + opt = pyo.SolverFactory("ipopt") + if max_iter is not None: + opt.options["max_iter"] = max_iter + opt.solve(m) + return m + + def test_sign_table_matches_ipopt(self): + for kind, expected in IPOPT_DUAL_AT_OPTIMUM.items(): + with self.subTest(kind=kind): + m = self._solve(_model(kind)) + self.assertAlmostEqual(pyo.value(m.dual[m.c]), expected, places=5) + + def test_tight_at_convergence(self): + for kind in IPOPT_DUAL_AT_OPTIMUM: + with self.subTest(kind=kind): + q = certified_lower_bound(self._solve(_model(kind)), eps_rel=0.0) + self.assertLessEqual(q, OPT + 1e-9) + self.assertLess(OPT - q, 1e-6) + + def test_valid_under_truncated_solves(self): + # Validity must hold at every truncation level; this is the property a + # lost minus sign would break by producing a bound *above* the optimum. + for max_iter in (1, 2, 3, 5, 8): + with self.subTest(max_iter=max_iter): + q = certified_lower_bound( + self._solve(_model("le"), max_iter=max_iter), eps_rel=0.0 + ) + if q is not None: + self.assertLessEqual(q, OPT + 1e-9) + + def test_converged_bound_is_at_least_as_tight_as_truncated(self): + # Deliberately not asserting pairwise monotonicity across max_iter: the + # iterate path is a solver detail and may differ by linear-solver build. + converged = certified_lower_bound(self._solve(_model("le")), eps_rel=0.0) + for max_iter in (1, 2, 3, 5): + with self.subTest(max_iter=max_iter): + q = certified_lower_bound( + self._solve(_model("le"), max_iter=max_iter), eps_rel=0.0 + ) + if q is not None: + self.assertLessEqual(q, converged + 1e-9) + + def test_solver_objective_value_can_be_an_invalid_bound(self): + # The motivating failure: at a badly truncated solve the returned + # objective value sits *above* the optimum, so "just use f(vhat)" is not + # an outer bound at all, while the certificate stays valid. + m = self._solve(_model("le"), max_iter=1) + self.assertGreater(pyo.value(m.obj), OPT) + q = certified_lower_bound(m, eps_rel=0.0) + self.assertTrue(q is None or q <= OPT + 1e-9) + + +# --------------------------------------------------------------------------- +# Ill-conditioning +# +# The question this answers is whether a solve that goes badly because the +# problem is badly conditioned can produce an *invalid* bound rather than +# merely a loose one. It cannot, and the reason is structural: the certificate +# is an evaluation, not a solve. A condition number measures how much a linear +# solve amplifies error, and `certified_lower_bound` never solves anything -- it +# evaluates phi and its gradient at whatever point ipopt stopped at. The +# tangent inequality it rests on is pointwise and exact for any convex phi at +# any vhat, with no error constant, so kappa has nowhere to enter. What +# conditioning does change is how far vhat lands from optimal and how large the +# gradient there is, and both of those show up in the correction term as +# looseness. +# +# The test problem is a Hilbert-matrix QP, the standard ill-conditioned test +# case: min 1/2 x'Hx - 1'x s.t. sum(x) <= 5, x in [-10,10]^n, with +# H_ij = 1/(i+j+1). It is convex (H is positive definite) so the hypotheses +# hold exactly, while cond(H) is about 1.6e13 at n=10 and 3.5e17 at n=16 -- +# past the point where a double-precision solve can converge properly. +# --------------------------------------------------------------------------- + +HILBERT_BOX = 10.0 +HILBERT_RHS = 5.0 + + +def _hilbert_qp(n=10, rowscale=1.0, offset=0.0): + """min 1/2 x'Hx - 1'x + offset s.t. rowscale*sum(x) <= rowscale*5. + + `rowscale` leaves the feasible set alone and rescales the row, which + rescales the multiplier ipopt reports by the same factor -- 1e-8 drives it + to ~1e8 and makes the `lam^T g` term in phi dominate an objective of order + one, which is the cancellation regime. `offset` does the same to the `f` + term. Both are ways of making the arithmetic hard without touching the + mathematics. + """ + m = pyo.ConcreteModel() + m.I = pyo.RangeSet(0, n - 1) + m.x = pyo.Var(m.I, bounds=(-HILBERT_BOX, HILBERT_BOX), initialize=0.0) + m.obj = pyo.Objective( + expr=0.5 * sum(m.x[i] * m.x[j] / (i + j + 1) for i in range(n) for j in range(n)) + - sum(m.x[i] for i in range(n)) + + offset + ) + m.c = pyo.Constraint( + expr=rowscale * sum(m.x[i] for i in range(n)) <= rowscale * HILBERT_RHS + ) + m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT) + return m + + +def _hilbert_objective(xs, offset): + n = len(xs) + return ( + 0.5 * sum(xs[i] * xs[j] / (i + j + 1) for i in range(n) for j in range(n)) + - sum(xs) + + offset + ) + + +def _rigorous_upper_bound(xs, offset): + """f at a point PROVABLY in the feasible set, hence >= the true optimum. + + Comparing the certificate against `f(vhat)` from a converged solve is the + trap this exists to avoid: ipopt's vhat carries up to `constr_viol_tol` of + infeasibility, so `f(vhat)` can sit *below* the optimum and manufacture an + apparent violation where there is none. Here a uniform downshift is + bisected until the shifted, clipped point satisfies the constraint exactly + as evaluated -- clipping keeps it in the box and the shifted sum is + monotone in the shift, so the result is feasible by construction, whatever + the solver did. + """ + def clip(t): + return [min(HILBERT_BOX, max(-HILBERT_BOX, v - t)) for v in xs] + + lo, hi = 0.0, 2.0 * HILBERT_BOX + abs(HILBERT_RHS) + 1.0 + for _ in range(200): + mid = 0.5 * (lo + hi) + if sum(clip(mid)) > HILBERT_RHS: + lo = mid + else: + hi = mid + feasible = clip(hi) + assert sum(feasible) <= HILBERT_RHS + assert max(abs(v) for v in feasible) <= HILBERT_BOX + return _hilbert_objective(feasible, offset) + + +@unittest.skipUnless(ipopt_available, "ipopt is not available") +class TestIllConditioning(unittest.TestCase): + """Validity must survive a solve that goes badly for numerical reasons.""" + + # (label, n, rowscale, offset) + VARIANTS = [ + ("cond 1.6e13", 10, 1.0, 0.0), + ("tiny multiplier", 10, 1e8, 0.0), + ("multiplier ~1e8", 10, 1e-8, 0.0), + ("cond 3.5e17", 16, 1.0, 0.0), + ("offset 1e12", 10, 1.0, 1e12), + ] + TRUNCATIONS = (1, 2, 3, 5, 10, None) + + @staticmethod + def _solve(m, max_iter=None): + opt = pyo.SolverFactory("ipopt") + if max_iter is not None: + opt.options["max_iter"] = max_iter + opt.solve(m) + return m + + def _reference(self, n, rowscale, offset): + m = self._solve(_hilbert_qp(n, rowscale, offset)) + return _rigorous_upper_bound([pyo.value(m.x[i]) for i in range(n)], offset) + + def test_bound_never_exceeds_a_feasible_objective(self): + for label, n, rowscale, offset in self.VARIANTS: + upper = self._reference(n, rowscale, offset) + # Both sides are computed in double precision, so the comparison + # gets a few-thousand-ulp relative allowance. That is the + # arithmetic caveat dual_certificate.py documents, not slack in the + # theorem: at an offset of 1e12 a double resolves about 1e-4 + # absolute, and an exact comparison there would be testing the + # floating-point unit rather than the certificate. + slop = 1e-12 * (1.0 + abs(upper)) + for max_iter in self.TRUNCATIONS: + with self.subTest(variant=label, max_iter=max_iter): + m = self._solve(_hilbert_qp(n, rowscale, offset), max_iter) + q = certified_lower_bound(m, eps_rel=0.0) + if q is None: + continue + self.assertTrue(math.isfinite(q)) + self.assertLessEqual(q, upper + slop) + + def test_severe_row_scaling_really_does_produce_a_huge_multiplier(self): + # Without this the cancellation variant above could silently stop + # exercising cancellation -- a passing test that tests nothing. + m = self._solve(_hilbert_qp(10, rowscale=1e-8)) + self.assertGreater(abs(pyo.value(m.dual[m.c])), 1e6) + + def test_conditioning_costs_tightness_not_validity(self): + # The well-scaled case closes to the last few digits; the badly scaled + # one does not close at all. Both stay valid, which is the whole + # claim: conditioning moves the bound down, never up. + good_n, good_scale = 10, 1.0 + good_upper = self._reference(good_n, good_scale, 0.0) + good = certified_lower_bound( + self._solve(_hilbert_qp(good_n, good_scale)), eps_rel=0.0) + self.assertIsNotNone(good) + self.assertLessEqual(good, good_upper + 1e-12 * (1.0 + abs(good_upper))) + self.assertLess(good_upper - good, 1e-6 * (1.0 + abs(good_upper))) + + bad_upper = self._reference(10, 1e-8, 0.0) + bad = certified_lower_bound( + self._solve(_hilbert_qp(10, rowscale=1e-8)), eps_rel=0.0) + self.assertIsNotNone(bad) + self.assertLessEqual(bad, bad_upper + 1e-12 * (1.0 + abs(bad_upper))) + + +if __name__ == "__main__": + unittest.main() diff --git a/mpisppy/tests/test_flexible_rank_cli.py b/mpisppy/tests/test_flexible_rank_cli.py index 65d501f16..8d352ab95 100644 --- a/mpisppy/tests/test_flexible_rank_cli.py +++ b/mpisppy/tests/test_flexible_rank_cli.py @@ -28,6 +28,7 @@ def _full_cfg(): cfg.fwph_args() cfg.lagrangian_args() cfg.gapper_args("lagrangian") + cfg.ipopt_outer_bound_args() cfg.ph_dual_args() cfg.relaxed_ph_args() cfg.ph_xfeas_spoke_args() @@ -57,6 +58,7 @@ def test_rank_ratio_args_default_to_one(self): self.assertEqual(cfg.ph_dual_rank_ratio, 1.0) self.assertEqual(cfg.relaxed_ph_rank_ratio, 1.0) self.assertEqual(cfg.subgradient_rank_ratio, 1.0) + self.assertEqual(cfg.ipopt_outer_bound_rank_ratio, 1.0) def test_build_spoke_list_injects_rank_ratio(self): cfg = _full_cfg() @@ -77,6 +79,8 @@ def test_build_spoke_list_injects_rank_ratio(self): cfg.relaxed_ph_rank_ratio = 5.0 cfg.subgradient = True cfg.subgradient_rank_ratio = 0.125 + cfg.ipopt_outer_bound = True + cfg.ipopt_outer_bound_rank_ratio = 6.0 scenario_creator = farmer.scenario_creator scenario_denouement = farmer.scenario_denouement @@ -92,7 +96,7 @@ def test_build_spoke_list_injects_rank_ratio(self): # exactly the enabled spokes, each carrying its requested ratio ratios = sorted(d["rank_ratio"] for d in spokes) self.assertEqual(ratios, sorted([0.5, 0.25, 2.0, 4.0, - 8.0, 3.0, 5.0, 0.125])) + 8.0, 3.0, 5.0, 0.125, 6.0])) # and every spoke dict got an explicit rank_ratio self.assertTrue(all("rank_ratio" in d for d in spokes)) diff --git a/mpisppy/tests/test_ipopt_outer_bound.py b/mpisppy/tests/test_ipopt_outer_bound.py new file mode 100644 index 000000000..88014292a --- /dev/null +++ b/mpisppy/tests/test_ipopt_outer_bound.py @@ -0,0 +1,370 @@ +############################################################################### +# mpi-sppy: MPI-based Stochastic Programming in PYthon +# +# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for +# Sustainable Energy, LLC, The Regents of the University of California, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for +# full copyright and license information. +############################################################################### +# Tests for the ipopt_outer_bound spoke. +# +# python -m pytest mpisppy/tests/test_ipopt_outer_bound.py +# mpiexec -np 2 python -m mpi4py -m pytest mpisppy/tests/test_ipopt_outer_bound.py +# +# The wiring tests need neither a solver nor MPI: option routing is the part +# most likely to break silently, and it is checkable by inspecting the spoke +# dict the factory builds. The end-to-end test needs both Ipopt and two ranks +# and skips cleanly without them. + +import math +import unittest + +import pyomo.environ as pyo + +import mpisppy.tests.examples.farmer as farmer +import mpisppy.utils.cfg_vanilla as vanilla +from mpisppy.utils import config +from mpisppy.spin_the_wheel import WheelSpinner +from mpisppy.utils.dual_certificate import CertificateError +from mpisppy.tests.utils import announce_hsl_if_used, get_solver + +from mpi4py import MPI + +comm = MPI.COMM_WORLD + +ipopt_available = pyo.SolverFactory("ipopt").available(exception_flag=False) + +if ipopt_available: + announce_hsl_if_used() +mip_available, mip_solver_name, *_ = get_solver() + + +def _reports_dual_bound(name): + """True if `name` actually fills in a dual bound on a solved LP. + + Being available is not enough, and neither is solving to optimality. cbc on + the CI runner returns `status=ok, TerminationCondition=optimal` and leaves + Problem[0].Lower_bound empty, so the Lagrangian spoke gets nothing to send + and its reported bound stays nan. Asking the solver directly is the only + honest test; anything else guesses. + """ + try: + if not pyo.SolverFactory(name).available(exception_flag=False): + return False + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 10), initialize=0.0) + m.o = pyo.Objective(expr=2.0 * m.x, sense=pyo.minimize) + m.c = pyo.Constraint(expr=m.x >= 1) + results = pyo.SolverFactory(name).solve(m, load_solutions=False) + bound = results.Problem[0].Lower_bound + return bound is not None and float(bound) == float(bound) # not NaN + except Exception: + return False + + +def _dual_bound_solver(): + """A solver that actually reports a dual bound, for the Lagrangian spoke. + + Ipopt reports none -- that is the entire reason this spoke exists -- so the + comparison in TestAgreesWithLagrangian needs a second solver. cbc is tried + as a fallback because it ships in the same idaes-ext bundle that supplies + Ipopt, but it is only used if it demonstrably reports a bound here. + """ + candidates = [] + if mip_available: + # The persistent interfaces need set_instance before a solve, which the + # probe below does not do, so try the plain name as well. + candidates += [mip_solver_name, mip_solver_name.replace("_persistent", "")] + # glpk and cbc are both plausible in CI: glpk is a small apt package and cbc + # ships in the idaes-ext bundle alongside Ipopt. Neither is assumed to work + # -- each is probed. (glpk cannot handle a PH proximal term, but the + # Lagrangian spoke runs with attach_prox=False, so it is fine here.) + candidates += ["glpk", "cbc"] + for name in candidates: + if name and _reports_dual_bound(name): + return name + return None + + +dual_bound_solver_name = _dual_bound_solver() + +# The three-scenario farmer optimum, from an EF solve. farmer is linear, hence +# convex, so it is inside this spoke's assumptions and its bound must not exceed +# this value. +FARMER_EF_OPT = -108390.0 + + +def _cfg(num_scens=3, hub_solver="ipopt"): + cfg = config.Config() + cfg.num_scens_required() + cfg.popular_args() + cfg.two_sided_args() + cfg.ph_args() + cfg.lagrangian_args() + cfg.ipopt_outer_bound_args() + cfg.num_scens = num_scens + cfg.max_iterations = 5 + cfg.default_rho = 1.0 + cfg.solver_name = hub_solver + return cfg + + +def _beans(cfg): + all_scenario_names = farmer.scenario_names_creator(cfg.num_scens) + kwargs = farmer.kw_creator(cfg) + beans = (cfg, farmer.scenario_creator, farmer.scenario_denouement, + all_scenario_names) + return beans, kwargs + + +def _spoke_options(cfg, **kw): + beans, kwargs = _beans(cfg) + spoke = vanilla.ipopt_outer_bound_spoke(*beans, + scenario_creator_kwargs=kwargs, **kw) + return spoke["opt_kwargs"]["options"] + + +class TestConfigSurface(unittest.TestCase): + + def test_flags_exist_with_expected_defaults(self): + cfg = _cfg() + self.assertFalse(cfg.ipopt_outer_bound) + self.assertEqual(cfg.ipopt_outer_bound_rank_ratio, 1.0) + self.assertEqual(cfg.ipopt_outer_bound_cushion, 1e-9) + # Scoped to Ipopt, but the name is still overridable. + self.assertIn("ipopt_outer_bound_solver_name", cfg) + + def test_no_mipgap_flags(self): + # Ipopt is not a branch-and-bound solver; offering mipgap flags would + # imply otherwise. + cfg = _cfg() + self.assertNotIn("ipopt_outer_bound_starting_mipgap", cfg) + self.assertNotIn("ipopt_outer_bound_iter0_mipgap", cfg) + + +class TestFactoryWiring(unittest.TestCase): + + def test_solver_defaults_to_ipopt(self): + # Even when the hub runs something else entirely, the spoke must land + # on ipopt rather than inheriting -- its own setup guard would reject + # anything else. + options = _spoke_options(_cfg(hub_solver="gurobi")) + self.assertEqual(options["solver_name"], "ipopt") + + def test_explicit_solver_name_is_honored(self): + cfg = _cfg() + cfg.ipopt_outer_bound_solver_name = "ipopt_v2" + self.assertEqual(_spoke_options(cfg)["solver_name"], "ipopt_v2") + + def test_cushion_is_threaded_through(self): + cfg = _cfg() + cfg.ipopt_outer_bound_cushion = 1e-7 + self.assertEqual( + _spoke_options(cfg)["ipopt_outer_bound_cushion"], 1e-7) + + def test_global_solver_options_do_not_leak(self): # noqa: D401 + # The point of this test: Ipopt hard-fails on an unrecognized keyword + # rather than ignoring it, so inheriting the global --solver-options + # (meant for the hub's MIP solver) would kill this spoke on its first + # solve, with an error naming Ipopt rather than the option routing. + cfg = _cfg() + cfg.solver_options = "mipgap=0.01" + options = _spoke_options(cfg) + self.assertNotIn("mipgap", options["iter0_solver_options"]) + self.assertNotIn("mipgap", options["iterk_solver_options"]) + self.assertEqual(options["solver_options_layers"], []) + + def test_other_spokes_still_inherit_global_options(self): + # The contrast that makes the previous test meaningful: not inheriting + # is special to this spoke, not a change in how spokes work. + cfg = _cfg() + cfg.solver_options = "mipgap=0.01" + beans, kwargs = _beans(cfg) + lag = vanilla.lagrangian_spoke(*beans, scenario_creator_kwargs=kwargs) + self.assertIn("mipgap", + lag["opt_kwargs"]["options"]["iter0_solver_options"]) + + def test_max_solver_threads_does_not_leak(self): + # --max-solver-threads is re-applied by apply_solver_specs *after* the + # factory clears the global layers, so clearing alone is not enough. + # Ipopt has no `threads` option and translate_solver_options has no + # mapping for it, so it would reach the solver verbatim and hard-fail + # the spoke's first solve -- taking the whole run with it. + cfg = _cfg() + cfg.max_solver_threads = 2 + options = _spoke_options(cfg) + self.assertNotIn("threads", options["iter0_solver_options"]) + self.assertNotIn("threads", options["iterk_solver_options"]) + for layer in options["solver_options_layers"]: + self.assertNotIn("threads", layer["options"]) + + def test_max_solver_threads_stripped_but_spoke_options_kept(self): + # Stripping the cap must not take the spoke's own options with it. + cfg = _cfg() + cfg.max_solver_threads = 2 + cfg.ipopt_outer_bound_solver_options = "max_iter=42" + options = _spoke_options(cfg) + self.assertNotIn("threads", options["iterk_solver_options"]) + self.assertEqual(options["iterk_solver_options"].get("max_iter"), 42) + + def test_per_spoke_solver_options_do_apply(self): + # Not inheriting the global layer must not mean ignoring the spoke's + # own options, which is how Ipopt settings are meant to arrive. + cfg = _cfg() + cfg.solver_options = "mipgap=0.01" + cfg.ipopt_outer_bound_solver_options = "max_iter=42" + options = _spoke_options(cfg) + self.assertEqual(options["iterk_solver_options"].get("max_iter"), 42) + self.assertNotIn("mipgap", options["iterk_solver_options"]) + + +class TestSetupGuards(unittest.TestCase): + """The guards that belong to the spoke rather than the certificate engine. + + Constructed without running a wheel: the guard reads self.opt.options, so a + lightweight stand-in exercises it without a solve. + """ + + def _guard_with_solver(self, solver_name): + from mpisppy.cylinders.ipopt_outer_bound import IpoptOuterBound + + class _Stub: + options = {"solver_name": solver_name} + local_scenarios = {} + + spoke = IpoptOuterBound.__new__(IpoptOuterBound) + spoke.opt = _Stub() + spoke.cylinder_rank = 0 + return spoke + + def test_non_ipopt_solver_is_rejected(self): + spoke = self._guard_with_solver("gurobi") + with self.assertRaisesRegex(CertificateError, "scoped to Ipopt"): + spoke._check_setup_guards() + + def test_ipopt_variants_are_accepted(self): + # ipopt_v2 and similar names still name Ipopt. + for name in ("ipopt", "ipopt_v2"): + with self.subTest(name=name): + self._guard_with_solver(name)._check_setup_guards() + + def test_missing_solver_name_is_rejected(self): + spoke = self._guard_with_solver(None) + with self.assertRaisesRegex(CertificateError, "scoped to Ipopt"): + spoke._check_setup_guards() + + +@unittest.skipUnless(ipopt_available, "ipopt is not available") +@unittest.skipUnless(comm.size == 2, "needs exactly two ranks") +class TestAgainstEFOptimum(unittest.TestCase): + """End-to-end: a PH hub on a MIP solver, this spoke on Ipopt. + + This also exercises the routing claim in the design -- the hub and the spoke + really do run different solvers on their own copies of the models -- and the + Ebound reduction across the spoke's rank. + """ + + def _spin(self, hub_solver): + cfg = _cfg(hub_solver=hub_solver) + beans, kwargs = _beans(cfg) + hub_dict = vanilla.ph_hub(*beans, scenario_creator_kwargs=kwargs) + spoke = vanilla.ipopt_outer_bound_spoke( + *beans, scenario_creator_kwargs=kwargs) + wheel = WheelSpinner(hub_dict, [spoke]) + wheel.spin() + return wheel + + def _assert_valid_and_useful(self, wheel): + if wheel.global_rank != 1: + return + bound = wheel.spcomm.bound + self.assertIsNotNone(bound) + # An outer bound on a minimization must not exceed the optimum. + self.assertLessEqual(bound, FARMER_EF_OPT + 1e-6) + # And it must be useful, not merely valid: farmer's Lagrangian bound + # sits in the same neighborhood as the optimum, so a wildly negative + # number would mean the certificate had collapsed. + self.assertGreater(bound, 2.0 * FARMER_EF_OPT) + + def test_bound_does_not_exceed_the_ef_optimum(self): + # farmer is an LP, so Ipopt can drive the hub too; this keeps the test + # runnable anywhere Ipopt is, with no MIP solver needed. + self._assert_valid_and_useful(self._spin("ipopt")) + + @unittest.skipUnless(mip_available, "no MIP solver available") + def test_hub_and_spoke_can_use_different_solvers(self): + # The routing claim in the design: each cylinder solves its own copy of + # the models with its own solver, and the only coupling is the numeric + # exchange of W and bounds. Here the hub runs a MIP solver while the + # spoke runs Ipopt. + self._assert_valid_and_useful(self._spin(mip_solver_name)) + + +@unittest.skipUnless(ipopt_available, "ipopt is not available") +@unittest.skipUnless(comm.size == 2, "needs exactly two ranks") +@unittest.skipUnless(dual_bound_solver_name, "no dual-bound-reporting solver") +class TestAgreesWithLagrangian(unittest.TestCase): + """The strongest correctness check available: on a linear problem, compare + this spoke's bound against the ordinary Lagrangian spoke's. + + farmer is an LP, so an LP solver's dual bound *is* the Lagrangian dual value + -- exact, and arrived at by a completely different route than the tangent + plane over the variable box that this spoke computes from Ipopt's duals. + Two independent computations of the same quantity is a much sharper test + than "the bound does not exceed the optimum", which a badly broken + certificate could still pass by being very negative. + + Both legs use the same hub, same rho and same iteration count, so the hub + walks the same W trajectory and the two spokes are asked for a bound on the + same relaxations. Spoke bounds do not feed back into W. + """ + + TOL = 1e-2 # measured agreement is ~1.2e-4 + + def _bound_from(self, spoke_factory, **cfg_overrides): + cfg = _cfg(hub_solver="ipopt") + for k, v in cfg_overrides.items(): + setattr(cfg, k, v) + beans, kwargs = _beans(cfg) + hub_dict = vanilla.ph_hub(*beans, scenario_creator_kwargs=kwargs) + spoke = spoke_factory(*beans, scenario_creator_kwargs=kwargs) + wheel = WheelSpinner(hub_dict, [spoke]) + wheel.spin() + return wheel + + def test_certificate_reproduces_the_lagrangian_bound(self): + lag = self._bound_from(vanilla.lagrangian_spoke, + lagrangian_solver_name=dual_bound_solver_name) + cert = self._bound_from(vanilla.ipopt_outer_bound_spoke) + + if lag.global_rank != 1: + return + lag_bound, cert_bound = lag.spcomm.bound, cert.spcomm.bound + self.assertIsNotNone(lag_bound) + self.assertIsNotNone(cert_bound) + # A spoke that never sent anything leaves its bound at nan. That means + # the comparison solver produced no dual bound after all, which is a + # broken premise for this test rather than a failure of the spoke under + # test -- say so instead of reporting a bogus mismatch. + if math.isnan(lag_bound): + self.skipTest( + f"{dual_bound_solver_name} reported no dual bound for the " + "Lagrangian spoke, so there is nothing to compare against") + self.assertFalse(math.isnan(cert_bound), + "ipopt_outer_bound sent no bound at all") + + # Two independent computations of the same number. + self.assertAlmostEqual(cert_bound, lag_bound, delta=self.TOL) + + # And in the safe direction: the certificate carries the box-correction + # term and the cushion, so it may be slightly looser than the exact LP + # dual bound but must never be more optimistic than it. + self.assertLessEqual(cert_bound, lag_bound + self.TOL) + + # Both remain valid outer bounds. + self.assertLessEqual(cert_bound, FARMER_EF_OPT + 1e-6) + self.assertLessEqual(lag_bound, FARMER_EF_OPT + 1e-6) + + +if __name__ == "__main__": + unittest.main() diff --git a/mpisppy/tests/utils.py b/mpisppy/tests/utils.py index 5e7b1c344..4776aa5dd 100644 --- a/mpisppy/tests/utils.py +++ b/mpisppy/tests/utils.py @@ -8,6 +8,9 @@ ############################################################################### +import os +import tempfile + import pyomo.environ as pyo from math import log10, floor @@ -50,3 +53,95 @@ def get_solver(persistent_OK=True): def round_pos_sig(x, sig=1): return round(x, sig-int(floor(log10(abs(x))))-1) + + +# --- HSL acknowledgement ----------------------------------------------------- +# +# Ipopt builds that link the Harwell Subroutine Library print, in their own +# banner, that "any publicity material resulting from use of the HSL codes +# within IPOPT must contain the acknowledgement: HSL, a collection of Fortran +# codes for large-scale scientific computation." Our test solves run with +# tee=False, so that banner never reaches the screen. These helpers put the +# acknowledgement back, and say which linear solver is actually in use -- +# worth knowing anyway, since the idaes-ext build defaults to ma27 rather than +# to MUMPS, and results can differ between the two. + +_HSL_ACK = ( + "HSL, a collection of Fortran codes for large-scale scientific " + "computation. See https://www.hsl.rl.ac.uk/" +) + +_hsl_probe_result = None # cache: (linear_solver_name, uses_hsl) +_hsl_announced = False + + +def ipopt_linear_solver(): + """Return (linear_solver_name, uses_hsl) for the ipopt on PATH. + + Ipopt names its linear solver in the banner it writes at the start of every + solve, so one trivial solve into a logfile is enough. Returns (None, False) + when ipopt is unavailable or the banner cannot be read. + """ + global _hsl_probe_result + if _hsl_probe_result is not None: + return _hsl_probe_result + + result = (None, False) + try: + if pyo.SolverFactory("ipopt").available(exception_flag=False): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(-10, 10), initialize=0.0) + m.o = pyo.Objective(expr=(m.x - 3) ** 2) + m.c = pyo.Constraint(expr=m.x <= 1) + fd, path = tempfile.mkstemp(suffix=".log") + os.close(fd) + try: + pyo.SolverFactory("ipopt").solve(m, logfile=path) + with open(path) as f: + text = f.read() + finally: + if os.path.exists(path): + os.remove(path) + name = None + for line in text.splitlines(): + if "running with linear solver" in line: + name = line.split("running with linear solver")[1] + name = name.strip().rstrip(".").split()[0] + break + result = (name, "compiled using HSL" in text) + except Exception: + # Never let a courtesy message break a test run. + result = (None, False) + + _hsl_probe_result = result + return result + + +def announce_hsl_if_used(): + """Print the HSL acknowledgement, once per run, if ipopt links HSL. + + Gated on MPI rank: these tests also run under mpiexec, and the project + convention is that such output comes from rank 0 only -- otherwise the + banner is emitted once per rank and interleaves with itself. + """ + global _hsl_announced + if _hsl_announced: + return + _hsl_announced = True + try: + from mpi4py import MPI + if MPI.COMM_WORLD.Get_rank() != 0: + return + except ImportError: + pass + name, uses_hsl = ipopt_linear_solver() + if not uses_hsl: + return + bar = "=" * 78 + print( + f"\n{bar}\n" + f"These tests solve with Ipopt built against HSL" + + (f" (linear solver: {name})" if name else "") + + f".\n{_HSL_ACK}\n{bar}", + flush=True, + ) diff --git a/mpisppy/utils/cfg_vanilla.py b/mpisppy/utils/cfg_vanilla.py index 9aa269a0d..6d7b35149 100644 --- a/mpisppy/utils/cfg_vanilla.py +++ b/mpisppy/utils/cfg_vanilla.py @@ -1172,6 +1172,75 @@ def lagrangian_spoke( return lagrangian_spoke +def ipopt_outer_bound_spoke( + cfg, + scenario_creator, + scenario_denouement, + all_scenario_names, + scenario_creator_kwargs=None, + rho_setter=None, + all_nodenames=None, + ph_extensions=None, + extension_kwargs=None, +): + from mpisppy.cylinders.ipopt_outer_bound import IpoptOuterBound + ipopt_ob_spoke = _PHBase_spoke_foundation( + IpoptOuterBound, + cfg, + scenario_creator, + scenario_denouement, + all_scenario_names, + scenario_creator_kwargs=scenario_creator_kwargs, + rho_setter=rho_setter, + all_nodenames=all_nodenames, + ph_extensions=ph_extensions, + extension_kwargs=extension_kwargs, + ) + + # This spoke does NOT inherit the global --solver-options layer, unlike + # every other spoke. Ipopt hard-fails on an unrecognized keyword rather + # than ignoring it, so an ordinary run -- a MIP solver and its options for + # the hub, this spoke attached alongside -- would kill the spoke on its + # first solve with an error naming Ipopt rather than the option routing. + # Ipopt-specific settings come in through --ipopt-outer-bound-solver-options. + # Filtering the global dict against a list of Ipopt-known keywords was + # considered and rejected: the list would have to track Ipopt releases, and + # silently dropping an option the user set is worse than never applying it. + options = ipopt_ob_spoke["opt_kwargs"]["options"] + options["iter0_solver_options"] = dict() + options["iterk_solver_options"] = dict() + options["solver_options_layers"] = [] + + apply_solver_specs("ipopt_outer_bound", ipopt_ob_spoke, cfg) + + # apply_solver_specs ends by re-applying --max-solver-threads as a + # system-level cap, *after* the reset above, so clearing the layers first is + # not enough. Ipopt has no `threads` option and translate_solver_options has + # no mapping for it, so it would reach the solver verbatim and hard-fail the + # spoke's first solve -- exactly the leak this factory is trying to prevent, + # reintroduced by a later step. Strip it here, where it is unambiguously + # wrong: the cap is meaningful for the MIP solvers it was written for. + for _key in ("iter0_solver_options", "iterk_solver_options"): + options[_key].pop("threads", None) + _stripped = [] + for _layer in options["solver_options_layers"]: + _layer["options"].pop("threads", None) + if _layer["options"]: + _stripped.append(_layer) + options["solver_options_layers"] = _stripped + + # apply_solver_specs only sets solver_name when the per-spoke flag was + # given; without this the spoke would silently inherit the global solver, + # which the setup guard then rejects. The spoke is Ipopt-scoped, so ipopt + # is the sensible default rather than something the user must repeat. + if not cfg.get("ipopt_outer_bound_solver_name"): + options["solver_name"] = "ipopt" + options["ipopt_outer_bound_cushion"] = cfg.ipopt_outer_bound_cushion + add_ph_tracking(ipopt_ob_spoke, cfg, spoke=True) + + return ipopt_ob_spoke + + def reduced_costs_spoke( cfg, scenario_creator, diff --git a/mpisppy/utils/config.py b/mpisppy/utils/config.py index 9b3f7a2e3..a49ce5d3a 100644 --- a/mpisppy/utils/config.py +++ b/mpisppy/utils/config.py @@ -1040,6 +1040,38 @@ def lagrangian_args(self): default=False) + def ipopt_outer_bound_args(self): + + self.add_to_config('ipopt_outer_bound', + description="have an ipopt_outer_bound spoke " + "(certified Lagrangian outer bound for " + "convex NLP subproblems; see spokes.rst)", + domain=bool, + default=False) + + self.add_to_config('ipopt_outer_bound_rank_ratio', + description="MPI ranks for the ipopt_outer_bound " + "spoke relative to the hub (flexible rank " + "assignments; default 1.0 = equal)", + domain=float, + default=1.0) + + # No add_mipgap_specs: Ipopt is not a branch-and-bound solver and has no + # mip gap. Offering the flags would suggest otherwise. + # The spoke is scoped to Ipopt, so its solver defaults to ipopt when + # this is left unset (applied in cfg_vanilla.ipopt_outer_bound_spoke; + # add_solver_specs itself defaults every solver name to None). + self.add_solver_specs("ipopt_outer_bound") + + self.add_to_config('ipopt_outer_bound_cushion', + description="relative cushion subtracted from the " + "certified bound: report q - eps*(1+|q|). " + "Last-bit hygiene against floating point, " + "not a proof-carrying margin; 0 disables", + domain=float, + default=1e-9) + + def reduced_costs_args(self): self.add_to_config('reduced_costs', diff --git a/mpisppy/utils/dual_certificate.py b/mpisppy/utils/dual_certificate.py new file mode 100644 index 000000000..3b69053be --- /dev/null +++ b/mpisppy/utils/dual_certificate.py @@ -0,0 +1,326 @@ +############################################################################### +# mpi-sppy: MPI-based Stochastic Programming in PYthon +# +# Copyright (c) 2024, Lawrence Livermore National Security, LLC, Alliance for +# Sustainable Energy, LLC, The Regents of the University of California, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md for +# full copyright and license information. +############################################################################### +# A *certified* lower bound for a convex NLP, computed from an already-solved +# Pyomo model and its constraint duals. +# +# The problem is +# +# L = min f(v) s.t. g(v) <= 0, h(v) == 0, v in B = [lo, hi] +# +# and the point of this module is to return a number that is <= L by a theorem, +# not by an assertion that the solver converged. A solver's returned objective +# value is the value at a point, hence >= L: an inner bound, the wrong +# direction. Lagrangian weak duality gives the right direction and needs no +# convergence assumption at all: +# +# phi(v) = f(v) + lam^T g(v) + mu^T h(v) +# q(lam, mu) = inf_{v in B} phi(v) <= L for ANY lam >= 0, ANY mu +# +# The remaining trap is that q is an *infimum*. Solving for it with an NLP +# solver returns a point, and the value there is >= q -- the wrong direction +# again. What closes it is convexity: phi is convex on B, so its tangent at any +# vhat in B lies below it, and minimizing that tangent over a box is separable +# and closed-form: +# +# q(lam, mu) >= phi(vhat) + sum_i min_{v_i in [lo_i, hi_i]} d_i phi * (v_i - vhat_i) +# +# One gradient evaluation and a loop over variables. No second solve, no +# tolerance argument, and vhat need not even be feasible. +# +# At an exact KKT point the correction term is identically zero (stationarity +# makes grad phi = z_L - z_U, which is zero on interior components and points +# the wrong way to help on active ones), so the bound collapses to f(vhat) and +# strong duality is recovered exactly. The correction is precisely the price of +# inexactness, and it measures itself. +# +# CONVEXITY IS LOAD-BEARING. If f or any component of g is non-convex over B, +# the tangent is not an underestimator and the returned number is simply wrong. +# +# Note carefully that the requirement is on the CANONICAL g, not on the +# constraint body as written, and the two differ by a sign on a `>=` row: +# +# body <= upper -> g = body - upper needs the body CONVEX +# body >= lower -> g = lower - body needs the body CONCAVE +# lo <= body <= up-> both of the above needs the body AFFINE +# body == rhs -> h = body - rhs needs the body AFFINE +# +# So `x**2 <= 4` is fine and `x**2 >= 1` is not, even though both are written +# with a convex body. This is easy to get backwards: `x**2 >= 1` looks like an +# ordinary convex constraint and its feasible set is not convex at all. +# +# check_model_is_certifiable() rejects what is mechanically checkable -- the +# affine cases, since polynomial degree is decidable -- but convexity of a +# general nonlinear body is not, so on one-sided nonlinear rows it is the +# caller's assertion. +# +# By contrast a wrong multiplier -- bad sign convention, stale dual, clipped +# value -- can only make the bound loose, never wrong, because weak duality +# holds for any lam >= 0 and any mu. In particular an inexact solve costs +# nothing but tightness: for a convex model there are no non-global local +# minima, so "the solver returned a sub-optimal answer" can only mean it stopped +# short of converging, and the theorem above never asked it to converge. +# +# THE ARITHMETIC IS A DIFFERENT MATTER, and it is the one place where an +# ill-conditioned model can hurt. phi(vhat) is evaluated in double precision as +# f + sum lam_i g_i + sum mu_j h_j, whose rounding error is on the order of +# u * (|f| + sum |lam_i g_i| + sum |mu_j h_j|) -- driven by the size of the +# TERMS. The eps_rel cushion below is proportional to |qhat| instead, i.e. to +# the size of the RESULT. On a well-scaled model those track each other; on one +# where cancellation is severe -- multipliers of 1e8 against an objective of +# order one, say -- they do not, and the cushion can be the smaller of the two. +# Hence the honest description of eps_rel as hygiene rather than proof: a user +# who needs margin on a badly conditioned model should raise it. +# +# The correction term carries the same kind of error, and it is the half that +# can push the answer UP: the correction is -sum_i |d_i phi| * d_i, so a +# gradient component computed slightly small in magnitude makes it slightly less +# negative. The clipping above gives no protection here -- that argument is +# about the multipliers, and a perturbed gradient is not one-directional the way +# a perturbed lam is. +# +# None of this is amplified by conditioning as such. A condition number +# measures how much a LINEAR SOLVE magnifies error, and nothing here solves +# anything: phi and its gradient are evaluated at vhat, each component is +# compared with zero, an endpoint is chosen, and the results are summed. The +# tangent inequality itself holds exactly at every vhat with no error constant. +# So an ill-conditioned subproblem gives a loose bound -- a worse vhat and a +# bigger gradient both enlarge the correction -- and not a wrong one. +# +# The caller is responsible for handing over a model that really is the +# relaxation it wants bounded. In particular a Progressive Hedging proximal +# term, or nonanticipative variables fixed by an extension, make the model +# something other than the Lagrangian relaxation, and this module cannot detect +# either. + +import math + +import pyomo.environ as pyo +from pyomo.core.expr.calculus.derivatives import differentiate, Modes +from pyomo.core.expr.visitor import identify_variables +from pyomo.contrib.fbbt.fbbt import fbbt + +__all__ = [ + "CertificateError", + "check_model_is_certifiable", + "unbounded_variables", + "certified_lower_bound", +] + + +class CertificateError(RuntimeError): + """The model violates an assumption the certified bound depends on.""" + + +# How a solver's Pyomo `dual` suffix maps onto canonical multipliers for +# g(v) <= 0 and h(v) == 0. These were measured against analytically known +# multipliers, not assumed: min (x-3)^2 with the constraint active has true +# multiplier 4, and ipopt reports d = -4 for `body <= upper`, d = +4 for +# `body >= lower`, and d = -4 for `body == rhs`. +_SIGN_CONVENTIONS = { + # g = body - upper (from `body <= upper`) + "ipopt": { + "lam_upper": lambda d: max(-d, 0.0), + # g = lower - body (from `body >= lower`) + "lam_lower": lambda d: max(d, 0.0), + # h = body - rhs (from `body == rhs`) + "mu": lambda d: -d, + }, +} + + +def _active_objective(model): + objs = list( + model.component_data_objects(pyo.Objective, active=True, descend_into=True) + ) + if len(objs) != 1: + raise CertificateError( + f"expected exactly one active Objective, found {len(objs)}" + ) + return objs[0] + + +def check_model_is_certifiable(model): + """Raise CertificateError if `model` violates an assumption the + certificate depends on. + + Checks only what is mechanically checkable. Convexity of the objective and + of one-sided nonlinear inequality bodies is the caller's assertion and is + *not* checked -- including the direction of it: a `<=` row needs a convex + body and a `>=` row needs a CONCAVE one, because the canonical g negates + the body on a `>=`. See the module docstring. + """ + obj = _active_objective(model) + if obj.sense != pyo.minimize: + raise CertificateError( + "certified_lower_bound is minimize-only; maximization is handled by " + "mirroring the model before calling here" + ) + + discrete = [ + v.name + for v in model.component_data_objects(pyo.Var, active=True, descend_into=True) + if not v.fixed and not v.is_continuous() + ] + if discrete: + raise CertificateError( + "a convexity claim is definitionally false with discrete variables; " + f"these are not continuous: {', '.join(sorted(discrete)[:10])}" + + (" ..." if len(discrete) > 10 else "") + ) + + # Both cases below need the body affine, and polynomial degree decides that, + # so these are real checks rather than assertions. A one-sided nonlinear row + # is left to the caller: whether its body is convex (for <=) or concave (for + # >=) is not decidable here. + nonlinear_eq = [] + nonlinear_ranged = [] + for con in model.component_data_objects( + pyo.Constraint, active=True, descend_into=True + ): + two_sided = (not con.equality) and con.has_lb() and con.has_ub() + if not (con.equality or two_sided): + continue + degree = con.body.polynomial_degree() + if degree is None or degree > 1: + (nonlinear_eq if con.equality else nonlinear_ranged).append(con.name) + if nonlinear_eq: + raise CertificateError( + "a nonlinear equality makes mu^T h non-convex for one sign of mu, " + "which breaks the underestimator; offending constraints: " + f"{', '.join(sorted(nonlinear_eq)[:10])}" + + (" ..." if len(nonlinear_eq) > 10 else "") + ) + if nonlinear_ranged: + raise CertificateError( + "a two-sided constraint splits into g = body - upper AND " + "g = lower - body, so its body would have to be both convex and " + "concave -- i.e. affine -- for the underestimator to hold on both " + "rows; offending constraints: " + f"{', '.join(sorted(nonlinear_ranged)[:10])}" + + (" ..." if len(nonlinear_ranged) > 10 else "") + ) + + +def unbounded_variables(model, do_fbbt=True): + """Names of variables that still lack a finite bound, after optionally + tightening with feasibility-based bounds tightening. + + `fbbt` mutates `model` in place, tightening variable bounds using the + constraints. That is sound here and makes the certificate tighter: it + shrinks the box `B` without removing any feasible point. + + A non-empty return does not mean no bound is possible -- an unbounded + variable only defeats the certificate if its gradient component in phi is + nonzero -- but it does mean `certified_lower_bound` may return None. + """ + if do_fbbt: + fbbt(model) + return sorted( + v.name + for v in model.component_data_objects(pyo.Var, active=True, descend_into=True) + if not v.fixed and (v.lb is None or v.ub is None) + ) + + +def _lagrangian_expression(model, conv): + """phi(v) = f(v) + lam^T g(v) + mu^T h(v), with the multipliers read off the + model's `dual` suffix and canonicalized.""" + if not hasattr(model, "dual"): + raise CertificateError( + "model has no `dual` Suffix; attach " + "pyo.Suffix(direction=pyo.Suffix.IMPORT) before solving" + ) + dual = model.dual + + terms = [_active_objective(model).expr] + for con in model.component_data_objects( + pyo.Constraint, active=True, descend_into=True + ): + if con not in dual: + raise CertificateError( + f"no dual available for constraint {con.name}; the solve did not " + "import one" + ) + d = float(pyo.value(dual[con])) + body = con.body + if con.equality: + mu = conv["mu"](d) + if mu != 0.0: + terms.append(mu * (body - pyo.value(con.upper))) + continue + # A ranged constraint carries one dual for two rows; splitting it and + # applying both rules lands the magnitude on the active side and zero + # on the other. + if con.has_ub(): + lam = conv["lam_upper"](d) + if lam != 0.0: + terms.append(lam * (body - pyo.value(con.upper))) + if con.has_lb(): + lam = conv["lam_lower"](d) + if lam != 0.0: + terms.append(lam * (pyo.value(con.lower) - body)) + return sum(terms) + + +def certified_lower_bound(model, sign_convention="ipopt", eps_rel=1e-9): + """A number guaranteed <= the model's optimal value, or None. + + `model` must already be solved, with its `dual` Suffix populated and its + variables holding the returned point. Neither optimality nor feasibility of + that point is required -- a truncated solve yields a valid but loose bound. + + Returns None when the box minimization is unbounded below, which happens + when a variable with an infinite bound has a nonzero gradient component in + phi, and also when the arithmetic produces a non-finite result -- a NaN or + an infinity arriving from a diverged solve. None means "no bound this + time", never "-inf". + + `eps_rel` shaves a relative cushion off the result. At the default 1e-9 + this is last-bit hygiene, not a proof-carrying margin; pass 0.0 to get the + theorem's quantity exactly. + """ + try: + conv = _SIGN_CONVENTIONS[sign_convention] + except KeyError: + raise CertificateError( + f"unknown sign convention {sign_convention!r}; known: " + f"{sorted(_SIGN_CONVENTIONS)}" + ) + + phi = _lagrangian_expression(model, conv) + vlist = list(identify_variables(phi, include_fixed=False)) + + correction = 0.0 + if vlist: + grad = differentiate(phi, wrt_list=vlist, mode=Modes.reverse_numeric) + for v, g in zip(vlist, grad): + g = float(g) + if g == 0.0: + continue + vhat = pyo.value(v) + # min over [lo, hi] of a linear term goes to whichever end the + # gradient points away from. + bound = v.lb if g > 0.0 else v.ub + if bound is None: + return None + correction += g * (bound - vhat) + + qhat = pyo.value(phi) + correction + + # A solve that diverged or failed numerically can leave NaN in the point or + # in the duals, and NaN propagates silently through everything above. It + # must not leave this function. Returning it would be safe today only by + # accident -- NaN loses every comparison, so the hub's `new > old` update + # test happens to reject it -- and +inf, which an infinite multiplier can + # also produce, would be an outright invalid bound if anything ever did + # accept it. Both are "no bound", which this module already has a word for. + if not math.isfinite(qhat): + return None + + return qhat - eps_rel * (1.0 + abs(qhat)) diff --git a/run_coverage.bash b/run_coverage.bash index cbac43991..505303a9c 100755 --- a/run_coverage.bash +++ b/run_coverage.bash @@ -213,7 +213,12 @@ run_phase "serial unit tests (serial)" \ mpisppy/tests/test_prox_approx.py \ mpisppy/tests/test_sep_rho.py \ mpisppy/tests/test_reduced_costs_fixer.py \ - mpisppy/tests/test_slammer.py + mpisppy/tests/test_slammer.py \ + mpisppy/tests/test_dual_certificate.py \ + mpisppy/tests/test_ipopt_outer_bound.py + +run_phase "test_ipopt_outer_bound (mpiexec -np 2)" \ + mpiexec -np 2 coverage run --rcfile="$PROJ_DIR/.coveragerc" -m mpi4py -m pytest mpisppy/tests/test_ipopt_outer_bound.py -v run_phase "test_conf_int_farmer (spawns mpiexec)" \ coverage run --rcfile=.coveragerc mpisppy/tests/test_conf_int_farmer.py