Introduction

POUNCE is a general interior-point method, implemented in pure Rust — one numerical backbone that now spans nonlinear, conic/quadratic, and polynomial global optimization rather than a single problem class. Its nonlinear-programming core began as a faithful port of the Ipopt filter line-search method — the algorithm, console output, and option semantics follow upstream Ipopt closely enough that anyone used to reading ipopt logs can drop in pounce without relearning where the numbers live — and it has since grown into a family of solvers sharing that backbone:

• Nonlinear programming — the filter line-search interior-point method (the Ipopt port) plus an active-set SQP path, for general smooth problems

  min  f(x)
  s.t. g_L <= g(x) <= g_U
       x_L <=   x  <= x_U
  

  where f and g are twice-continuously-differentiable.

• Conic & quadratic — LP, convex QP, second-order (SOCP), positive-semidefinite (SDP), and the non-symmetric exponential and power cones, each solved to the global optimum.

• Global optimization — certified global optima for nonconvex polynomial problems via SOS / Lasserre relaxations. Nonconvex problems that are not polynomial get a local answer; POUNCE has no spatial branch-and-bound solver.

See Choosing a Solver for which solver fits which problem.

Pure Rust by default

The default build is pure Rust — no Fortran, no commercial solver, no system BLAS required. The bundled FERAL backend provides a sparse symmetric LDLᵀ factorization. The HSL MA57 backend is available behind the optional ma57 feature for users who have a license for libcoinhsl and have it installed (see Installation).

Status

Production-ready for the core IPM workflow. The algorithm-side core, NLP interface, line search, filter, barrier update (monotone + Mehrotra adaptive), KKT solve, restoration phase, AMPL .nl reader, the C ABI (pounce-cinterface), the Python wrapper (pounce-solver), and the CLI all solve a wide range of NLPs from the standard test suites (Hock-Schittkowski, CUTEst, Mittelmann ampl-nlp, CHO parameter estimation, gas/water network design). Sensitivity analysis (sIPOPT port), reduced-Hessian computation, the auxiliary-equality + FBBT presolve, and the active-set SQP path are all wired in and available behind option keys. Existing PyIpopt / cyipopt / JuMP / AMPL clients link against libpounce_cinterface in place of libipopt unchanged.

The conic and global solvers are wired end-to-end alongside the NLP core: the convex interior-point solver (pounce-convex) handles LP / QP, SOCP, exponential / power cones, and small SDPs — with a Conic Benchmark Format (.cbf) reader cross-checked against the CBLIB tier — and adds SOS / Lasserre polynomial global optimization (sos_minimize). These are reachable from the CLI, the Python package, and the JSON solve report. There is no spatial branch-and-bound solver for general factorable nonconvex problems — outside the polynomial case, nonconvex models are solved locally.

License

EPL-2.0, the same license as upstream Ipopt.

Where to go next

• Installation — build and install POUNCE.
• Quick Start — solve your first problem.
• Running Solves — the command-line driver in depth.
• Acknowledgments — the papers behind the algorithm.

Installation

Three routes, in the order most people want them:

With pip

pip install pounce-solver

Prebuilt wheels for Linux, macOS, and Windows (CPython 3.9+). No Rust toolchain is involved. This installs both interfaces at once:

pounce problem.nl                       # the CLI
python -c "import pounce; print(pounce.__version__)"

For Pyomo models:

pip install pyomo-pounce

import pyomo.environ as pyo
results = pyo.SolverFactory("pounce").solve(model)

Optional extras — none needed for a normal solve:

pip install "pounce-solver[jax]"     # pounce.jax autodiff frontend
pip install "pounce-solver[torch]"   # pounce.torch autodiff frontend
pip install "pounce-solver[viz]"     # debugger plots (pounce-dbg-viz)
pip install "pounce-solver[gams]"    # GAMS solver link — see gams.md

If the CLI will not start (GLIBC_2.39 not found)

Releases up to and including 0.9.0 bundled a Linux CLI built against a newer glibc than the wheel advertised, so pounce fails to exec on older distributions (Debian 12, Ubuntu 22.04, RHEL/Rocky/Alma 8 and 9, and most HPC images) with:

pounce/bin/pounce: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.39' not found

import pounce still works; only the CLI (and therefore Pyomo, which shells out to it) is affected. This is fixed from 0.10.0 on — the CLI is now built inside the manylinux container, and the published 0.10.0 wheel’s binary references nothing above GLIBC_2.16, under the manylinux2014 floor the wheel advertises. scripts/check-cli-portability.sh asserts that on every build. If you are pinned to 0.9.0 or earlier and hit this, upgrade, or use the container or a source build below.

With a container

No toolchain and nothing installed on the host:

docker run --rm -v "$PWD:/work" ghcr.io/jkitchin/pounce:latest problem.nl
apptainer pull pounce.sif docker://ghcr.io/jkitchin/pounce:0.11.0

Both images carry the CLI, the Python API, and the Pyomo plugin. See Docker & Containers for tags, bind mounts, and a Slurm example.

From source

Prerequisites

A stable Rust toolchain. Nothing else is needed for the default pure-Rust build. Install Rust via rustup:

curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh
source "$HOME/.cargo/env"

Verify the install:

rustc --version && cargo --version

Build

From the repository root:

make            # release build of the workspace
make test       # run all tests
make clippy     # lint
make doc        # rustdoc for the Rust API

Install

make install                          # installs to $HOME/.local
sudo make install PREFIX=/usr/local   # or system-wide

This drops the pounce binary into $PREFIX/bin and the libpounce_cinterface shared library into $PREFIX/lib. Make sure $HOME/.local/bin is on your PATH, then verify:

pounce --version

HSL MA57 backend (optional)

The default FERAL backend needs no external libraries. To build with the HSL MA57 linear solver instead, you need a CoinHSL install whose lib/ directory holds libcoinhsl. Point the COINHSL_DIR environment variable at it and build with the ma57 feature:

export COINHSL_DIR=/path/to/CoinHSL
cargo build -p pounce-cli --release --features ma57

The feature makes MA57 available; selecting it is a separate step, because linear_solver defaults to feral in every build:

pounce problem.nl linear_solver=ma57

Build CoinHSL from https://www.hsl.rl.ac.uk/ipopt/. MA57 is primarily useful for benchmarking against upstream Ipopt; the FERAL backend is the supported default for everyday use, and a build without --features ma57 never touches COINHSL_DIR.

The build embeds $COINHSL_DIR/lib as an rpath, so the resulting binary finds libcoinhsl and its own dependencies (openblas, metis, libgfortran, libgomp) without LD_LIBRARY_PATH or DYLD_LIBRARY_PATH. If you relocate the CoinHSL install afterwards, rebuild — or, on macOS, rewrite the path with install_name_tool -rpath <old> <new>, which works because the link reserves header padding for it.

The ma57_* options (ma57_pivtol, ma57_pivot_order, ma57_pre_alloc, and the rest — see Options) are honoured by this build, and can be scoped to the restoration sub-solve with a resto. prefix, e.g. resto.ma57_pivtol=0.5. Before issue #825 they were accepted and silently discarded.

One of them is a POUNCE addition rather than an Ipopt option: ma57_batched_backsolve lets the limited-memory correction hand MA57 several right-hand sides at once. It is off by default because turning it on perturbs the iterate by about one bit and therefore moves the trajectory — read that section before using it.

Using POUNCE as a Rust library

The workspace is a set of library crates (see Algorithm & Workspace for the layout). To browse the Rust API, build and open the rustdoc:

make doc        # generates target/doc

Docker & Containers

Prebuilt images carry the whole POUNCE surface — the pounce CLI, the Python API, and the Pyomo plugin — with no Rust toolchain and no pip install on your side. This is the path of least resistance on a cluster, where you often cannot install a toolchain anyway, and the quickest way to try POUNCE without touching your environment.

docker pull ghcr.io/jkitchin/pounce:latest
docker run --rm ghcr.io/jkitchin/pounce:latest --version

Which tag

Pin X.Y.Z for anything you intend to re-run months later — latest and X.Y both move under you.

Published images are cut from releases only. There are edge and sha-<short> tags in the registry from when every commit to main was built, but they are frozen at whatever commit last published them — do not read edge as “tip of main”. To run an unreleased fix, build the image yourself from a checkout; it takes minutes and needs no Rust toolchain on your side beyond Docker:

git clone https://github.com/jkitchin/pounce && cd pounce
make docker          # -> pounce:dev, compiled from the tree you checked out
docker run --rm pounce:dev --version

Running solves

The entrypoint is the pounce CLI, so arguments after the image name go straight to it:

docker run --rm ghcr.io/jkitchin/pounce:latest --list-problems
docker run --rm ghcr.io/jkitchin/pounce:latest --problem rosenbrock

Your own problems live on the host, so mount a directory. The working directory inside the image is /work:

docker run --rm -v "$PWD:/work" ghcr.io/jkitchin/pounce:latest \
  problem.nl print_level=5 tol=1e-10

That writes problem.sol next to problem.nl in the mounted directory — see Running Solves for the full option surface. The image runs as UID 1000 rather than root, so files it creates in a bind mount are not root-owned on the host. If your host UID differs, add --user "$(id -u):$(id -g)".

Python and Pyomo

Override the entrypoint to get a shell or an interpreter:

docker run --rm -it --entrypoint python ghcr.io/jkitchin/pounce:latest
docker run --rm -it --entrypoint bash -v "$PWD:/work" ghcr.io/jkitchin/pounce:latest

Both import pounce (with numpy and scipy) and Pyomo’s SolverFactory('pounce') work out of the box:

docker run --rm -v "$PWD:/work" --entrypoint python \
  ghcr.io/jkitchin/pounce:latest my_model.py

The optional extras are not installed — no JAX, no PyTorch, no plotly, no GAMS bindings. Add what you need in a derived image:

FROM ghcr.io/jkitchin/pounce:0.11.0
USER root
RUN pip install --no-cache-dir "pounce-solver[jax]==0.11.0"
USER pounce

On a cluster (Apptainer / Singularity)

Most HPC sites run Apptainer (formerly Singularity) rather than Docker, because it needs no daemon and no root. It pulls Docker images directly:

apptainer pull pounce.sif docker://ghcr.io/jkitchin/pounce:0.11.0
apptainer run pounce.sif problem.nl

Two differences from docker run are worth knowing before you write a job script:

• You are yourself. Apptainer runs the container as your own user, not the image’s, so output files land with your ownership and no --user flag is needed.
• $HOME and $PWD are already there. Apptainer bind-mounts them by default, so a .nl file in your submit directory is usually visible with no -B flag at all. Add -B /scratch:/scratch (or your site’s equivalent) for anything outside them.

Build the .sif once on a login node and reuse it — pulling on every array task hammers the registry and will get rate-limited. In a Slurm script:

#!/bin/bash
#SBATCH --job-name=pounce
#SBATCH --cpus-per-task=4

apptainer exec $HOME/images/pounce.sif \
  pounce $SLURM_SUBMIT_DIR/problem.nl --json-output result.json

Pin the digest rather than the tag if the run has to be reproducible years later — a tag can be re-pushed, a digest cannot:

apptainer pull pounce.sif docker://ghcr.io/jkitchin/pounce@sha256:<digest>

Building your own

You do not need the images to be published — both Dockerfiles are in the repository and build from a clone. From the repository root:

make docker          # compiles the current working tree -> pounce:dev
make docker-release  # installs the released wheels -> pounce:<version>

make docker is the one to reach for when testing a branch: it compiles whatever is checked out, and stamps the commit into pounce --about so the image can say what it contains. make docker-release needs no Rust toolchain and takes seconds. See docker/README.md for build arguments and the .dockerignore caveat.

What is not in the image

The HSL MA57 backend is absent, and cannot be added by us: CoinHSL is license-restricted and not redistributable. Passing linear_solver=ma57 will not work in a container. The pure-Rust FERAL backend is the default everywhere and needs no external libraries — see Installation if you have a CoinHSL license and want a local build with it.

The GAMS link is likewise absent, since it needs your own GAMS install and license on the host. See GAMS.

Quick Start

This page assumes POUNCE is built and on your PATH (see Installation).

Solve an AMPL .nl file

pounce problem.nl

This solves the problem and writes a sibling problem.sol next to the input, following the AMPL solver convention. The console output mirrors upstream ipopt’s banner, per-iteration table, and final summary.

Append KEY=VALUE pairs to override options — the syntax and semantics match the upstream Ipopt CLI:

pounce problem.nl print_level=8 max_iter=500 tol=1e-10

See Solver Options for details.

Try a built-in problem

POUNCE ships several self-contained test problems that exercise the full pipeline without parsing a .nl file (run pounce --list-problems for the full set):

pounce --list-problems
pounce --problem rosenbrock
pounce --problem quadratic

From Python

import numpy as np
from pounce import minimize

res = minimize(lambda x: ((x - 1) ** 2).sum(), x0=np.zeros(3))
print(res.fun, res.x)

See the Python API chapter for the full cyipopt-compatible interface.

From Pyomo

import pyomo_pounce  # registers 'pounce'
from pyomo.environ import SolverFactory

SolverFactory('pounce').solve(model)

See the Pyomo chapter for details.

Full help

pounce --help

Choosing a Solver

POUNCE is not a single solver but a small family of them sharing one numerical backbone. This page is the map: what each solver is, when to reach for it, and how they fit together.

The one-sentence version: convex and conic problems are solved to the global optimum; nonconvex problems are solved locally by default, or to a certified global optimum via the SOS (polynomial) and spatial branch-and-bound (general) paths. Every solver, whatever its flavor, ultimately factorizes a symmetric KKT system through the shared pounce-linsol layer, which in turn drives a pluggable backend (FERAL by default, HSL MA57 optionally).

The solvers at a glance

When to reach for the active-set QP. auto never selects it: a cold, one-shot convex QP goes to the interior-point path, which is materially more robust on that workload (137 of the 138 Maros-Mészáros problems, against substantially fewer for a cold active-set solve). That is the character of the method rather than a defect — an active-set iteration count is combinatorial in the size of the active set, while an interior-point count is nearly independent of problem size. Choose solver_selection=qp-active-set when you want an exact vertex solution, or when you are solving a sequence of similar QPs — MPC steps, branch-and-bound nodes, continuation — where the working set carries across solves and solve_parametric can trace the homotopy from the previous solution instead of starting over.

POUNCE has no spatial branch-and-bound solver for general factorable nonconvex NLPs — no solver_selection=global CLI route, no minimize_global Python entry point. The only certified-global path for a nonconvex problem is SOS / Lasserre, and it covers polynomials only.

When to choose each

General nonlinear program (the common case) → NLP filter-IPM

If your model has nonlinear objective or constraints and you don’t know (or can’t assume) convexity, this is the default and the most mature path. It is POUNCE’s port of Ipopt’s filter line-search interior-point method: robust on nonconvex problems, with a feasibility restoration phase for hard starts and exact or limited-memory Hessians. It returns a local KKT point — for a nonconvex problem there is no global guarantee.

• CLI: pounce model.nl (or a built-in problem).
• Python: the cyipopt-style Problem class, or the scipy-style minimize facade.
• Reach for limited-memory Hessians (hessian_approximation=limited-memory) when second derivatives are unavailable or expensive.

A sequence of related NLPs, or a stable active set → NLP active-set SQP

Selected with algorithm=active-set-sqp. It solves the NLP as a sequence of quadratic subproblems (handed to pounce-qp), which warm-starts extremely well when the active set is stable across solves — e.g. a parametric sweep or a control loop. For a single cold solve of a general NLP, prefer the filter-IPM.

Linear or convex quadratic program → Convex IPM (LP/QP)

If P ⪰ 0 (or P = 0 for an LP), use the convex interior-point solver: it returns the global optimum, detects primal/dual infeasibility, and offers warm-starting, batched and multiple-RHS solving, a build-once / solve-many QpFactorization handle, and post-optimal sensitivity (QpSensitivity — the sIPOPT analog). The CLI’s auto routing classifies an .nl and sends LP/convex-QP problems here automatically, and a .nl carrying the sIPOPT sens_* suffixes is now answered here rather than rerouted to the general NLP engine — see LP/QP routing for what still routes away and why.

• Python: pounce.qp.solve_qp (and solve_qp_batch, solve_qp_multi_rhs).

Second-order, exponential, or power cones → Convex IPM (conic)

The same convex solver handles conic programs: second-order cones, the exponential and power cones that express geometric programming, entropy / log-sum-exp, logistic models, and p-norm constraints, and the positive-semidefinite cone for small dense SDPs. Also global. This is the path to use when you can cast a nominally-nonconvex problem into a convex cone — you trade modeling effort for a global guarantee. (The PSD cone is self-scaled and runs on the symmetric driver; the exp/power cones run on the non-symmetric HSDE driver, so the two families can’t yet be mixed in one problem.)

A common special case routes here automatically: a convex quadratically-constrained QP (QCQP). When auto routing finds a convex-quadratic inequality ½xᵀHx + aᵀx + b ≤ 0 (H ⪰ 0), it reformulates each such constraint to one second-order cone (H = FᵀF) and sends the whole problem to the conic solver — no .cbf and no manual cone bookkeeping needed. This works from a .nl/Pyomo model on the CLI and from minimize() in Python (which probes each constraint’s Hessian and only routes when it can prove the feasible set is convex). See LP / QP Solver Routing.

• Python: pounce.qp.solve_socp(..., cones=[("exp", 3), ("pow", 0.5), ...]) for an explicit cone program, or just minimize(...) for a convex QCQP.
• CLI: a Conic Benchmark Format file, pounce model.cbf (see the CBLIB benchmark tier), or any convex-QCQP .nl under auto routing.

Nonconvex problem, global optimum required → SOS (polynomials only)

When the problem is genuinely nonconvex and a local optimum is not good enough, the one path to a certified global optimum is for polynomials:

• Polynomial objective/constraints → SOS / Lasserre (sos_minimize, or pounce.sos_minimize). A single semidefinite program certifies the global minimum (the largest γ with p − γ in the Putinar cone), and the global minimizers are recovered from the moment matrix — even multiple ones, via a facial-reduction step. Best for modest degree and dimension; the SDP grows with the relaxation order.

If the problem is nonconvex and not polynomial (exp/ln/trig), POUNCE cannot certify a global optimum. Reformulate into the convex cone library if you can; otherwise multistart the local NLP solver and accept that the result is uncertified.

See Global Optimization for the SOS path in depth, and for the multistart fallback.

Indefinite QP, or a QP inner-solver → Active-set QP

pounce-qp is a sparse parametric active-set solver that accepts an indefinite Hessian (via inertia control), with two-sided bounds and factorization-reuse across a homotopy. It is the engine behind the active-set SQP path, and is the right choice for MPC-style problems or any setting where you re-solve a slowly-changing QP many times. Use the convex IPM instead when P ⪰ 0 and you want a single robust solve with infeasibility certificates.

It is the only entry point here that will take an indefinite P; both the convex IPM and the conic solver refuse one, because without a PSD Hessian the IPM’s optimality test accepts a saddle point and reports it as optimal (issue #112). Ask for the engine by name:

from pounce.qp import solve_qp
r = solve_qp(P=P_indefinite, c=c, lb=lb, ub=ub, method="active-set")

pounce model.nl solver_selection=qp-active-set   # class: nonconvex QP

Two things to hold onto about the answer. It is local: a nonconvex QP can have many KKT points, and what you get is the one the active set walked to, the same guarantee minimize gives on a nonconvex NLP. It is at least a minimum, though, which took a fix — inertia control leaves the first-order conditions satisfied at a saddle, and until issue #848 the engine reported one as optimal (min ½xᵀ[[1,5],[5,1]]x over [−1,1]² came back at objective 0 against a true minimum of −4). From 0.11.0 two guards stand behind an optimal here. The engine tests the reduced Hessian on its working set’s null space, and when it finds a feasible direction of negative curvature it follows it — to a better point, to an unbounded verdict if nothing blocks it, or to an honest non-optimal status if it runs out of budget. The driver then screens the result by exhibition, refusing a verdict only where it can walk a direction and hand you a strictly better feasible point; that one is not confined to the working set’s null space, so it reaches negative curvature hidden behind a bound whose multiplier is exactly zero, which the first cannot see. What you do not get is a global minimum, or a guarantee in the case where neither guard concludes — the first-order verdict stands there, as it always did. And the constraints must be linear — the curvature this engine controls is the objective’s, so an indefinite objective over any quadratic row, convex row or not, is a nonconvex QCQP and goes to the NLP filter-IPM instead (see LP/QP routing).

How to override the automatic routing

The CLI classifies each .nl problem and picks a solver, but you can force the choice:

pounce model.nl solver_selection=auto          # default: classify, then route
pounce model.nl solver_selection=nlp           # filter-IPM (or active-set-sqp via algorithm=)
pounce model.nl solver_selection=lp-ipm        # convex LP interior-point
pounce model.nl solver_selection=qp-ipm        # convex QP interior-point
pounce model.nl solver_selection=socp          # conic interior-point (convex QCQP)
pounce model.nl solver_selection=qp-active-set # active-set QP

solver_selection is an ordinary POUNCE option, not a command-line flag: it is passed as a trailing KEY=VALUE pair (the ipopt CLI convention), and so also works from an options file, the pounce_options environment variable, or Pyomo’s solver.options. Forcing a value the problem class does not support is rejected with a message rather than silently ignored.

See LP / QP Solver Routing for how classification works and when it falls back to the more general solver.

The shared backbone

Every interior-point and active-set solver above assembles a symmetric KKT system and factorizes it through pounce-linsol. That trait layer is backend-agnostic:

• FERAL (pounce-feral) — a pure-Rust sparse symmetric LDLᵀ factorization. The default; no external dependencies.
• HSL MA57 (pounce-hsl) — the well-known Harwell solver via libcoinhsl, enabled with the ma57 build feature for large or ill-conditioned systems.

Because the backend is pluggable, the same solver code runs on either without change.

Cross-cutting layers

These are not solvers you select, but stages and tools the solvers share:

• Presolve (pounce-presolve) — an optional front-end that tightens bounds (feasibility-based bound tightening), removes redundant rows, and repairs LICQ degeneracies before the solve.
• Restoration (pounce-restoration) — the feasibility-recovery phase the filter-IPM enters when a step cannot reduce both infeasibility and the objective; pounce-l1penalty offers an ℓ₁-exact penalty reformulation for degenerate / LICQ-violating problems.
• Sensitivity — pounce-sensitivity gives sIPOPT-style parametric steps and reduced Hessians for the NLP. QpSensitivity does the same kind of thing for the convex QP and conic arms — the parametric step, fix-relax, path following, activity classification, and a face decomposition for every cone family — over a shared core (pounce-sens-core), so the two cannot drift on what a kink is. They are not the same surface: the NLP arm additionally has the corrector, the covariance/identifiability statistics, and the directional decision at a kink, and the two reduced Hessians are different computations behind one word (sIPOPT’s Schur route vs a null-space projection). See Sensitivity Analysis and The convex/conic solver.
• Cone library (pounce-convex) — nonnegative, second-order, exponential, power, and (for small dense problems) positive-semidefinite cones, so small SDPs solve as a convex class. The PSD cone cannot yet be mixed with the exponential/power cones in one problem (they use different drivers).
• Solve report — every path can emit the machine-readable pounce.solve-report/v1 JSON (status, iterations, residuals, timing). See JSON Solve Report.

Global vs. local — the honest summary

POUNCE settles a problem globally along two routes, and locally along one:

• Global by convexity — LP, convex QP, SOCP, and the exponential / power / PSD cone classes. Local is global, so a convex or conic reformulation buys the guarantee outright.
• Global by certificate (polynomials) — the SOS / Lasserre optimizer certifies the global minimum of a nonconvex polynomial from a single SDP; see Global Optimization.
• Local for general NLP — the filter-IPM and SQP paths converge to a KKT point, which for a nonconvex problem carries no global guarantee.

There is no third route: a nonconvex, non-polynomial problem gets a local answer, and no certificate.

Two practical levers for a “global” answer: modeling (cast as much as you can into the convex cone library) and, when that is not possible, the SOS / Lasserre optimizer for polynomials.

Running Solves

The pounce command-line driver solves built-in TNLPs and AMPL .nl files. Its console output mirrors upstream ipopt’s banner, per-iteration table, and final summary, so anyone used to reading ipopt logs can read pounce logs unchanged.

Basic usage

pounce problem.nl
pounce problem.nl print_level=8 max_iter=500 tol=1e-10
pounce problem.nl linear_solver=ma57            # with --features ma57
pounce problem.nl --options-file tuned.opt      # upstream-format options file
pounce problem.nl option_file_name=tuned.opt   # the Ipopt spelling; same thing
pounce problem.nl --no-options-file            # ignore ./pounce.opt, ./ipopt.opt

Trailing KEY=VALUE pairs follow the same syntax and semantics as the upstream Ipopt CLI; they override values loaded from the options file. With no options file named, ./pounce.opt or ./ipopt.opt is read if present, as ipopt reads ./ipopt.opt. See Solver Options.

Built-in problems

pounce --list-problems
pounce --problem quadratic
pounce --problem rosenbrock

• quadratic — min (x[0]-3)² + (x[1]-4)² (unconstrained, optimum (3, 4)).
• rosenbrock — min 100·(x[1]-x[0]²)² + (1-x[0])² (unconstrained, optimum (1, 1)).
• bounded-quadratic — quadratic with box bounds 0 ≤ x ≤ 2 (optimum at the upper corner (2, 2)).
• eq-quadratic — min x[0]² + x[1]² s.t. x[0] + x[1] = 1 (a single equality).
• circle — min x[0] s.t. x[0]² + x[1]² = 1 (a nonlinear equality).
• infeasible-eq — two contradictory equalities (x[0]+x[1]=1 and =2); exercises the infeasibility-detection path.

Run pounce --list-problems for the authoritative list.

Built-in problems have no .nl stub, so they only write a .sol file when --sol-output is given explicitly.

Degenerate / MPCC NLPs — the ℓ₁-exact penalty-barrier wrapper

For problems where the standard IPM thrashes in restoration because LICQ fails at the iterate (degenerate equalities, MPCC-like complementarity), enable the Thierry–Biegler ℓ₁-exact penalty-barrier wrapper:

pounce problem.nl l1_exact_penalty_barrier=yes

The wrapper turns every equality row c_i(x) = g_i into a slack-relaxed c_i(x) − p_i + n_i = g_i with (p_i, n_i) ≥ 0, augments the objective by ρ · Σ(p + n), and runs a Byrd–Nocedal–Waltz outer loop that escalates ρ until the slacks collapse (constraints satisfied) or saturate (locally infeasible problem detected). The user-visible (x*, λ*) are reported in the original variable space.

For everyday use, the simpler form is an auto-fallback:

pounce problem.nl l1_fallback_on_restoration_failure=yes

POUNCE first runs the standard solve. If it terminates in Restoration_Failed, Infeasible_Problem_Detected, Solved_To_Acceptable_Level, Maximum_Iterations_Exceeded, or Not_Enough_Degrees_Of_Freedom, the wrapper is invoked transparently and the result is promoted to Solve_Succeeded only if the retry succeeds. Otherwise the original status is preserved.

The tuning knobs are listed under Solver Options.

AMPL / Pyomo solver mode

AMPL drivers — and Pyomo’s ASL interface — invoke a solver as solver problem.nl -AMPL. Pass -AMPL to run pounce that way:

pounce problem.nl -AMPL

It changes nothing about the solve itself; it switches the process to the AMPL exit-code contract (see below), so the driver reads the termination from the .sol file rather than the exit status. The pyomo-pounce package builds on top of this mode.

Dual sign conventions

Two different quantities are in play, and they differ by a sign. Getting them confused is easy, so pounce names them differently:

They are related by marginal = −λ. The .sol writer performs this negation, so a .sol (and therefore Pyomo’s model.dual) carries shadow prices with the same sign Ipopt, glpk, CBC and CONOPT report. For min x₀²+x₁² s.t. x₀+x₁ = b the optimum is b²/2, so a .sol written at b = 2 reports +2.

The Python API’s mult_g keeps the Lagrange-multiplier convention, which matches cyipopt and satisfies the stationarity identity above directly. If you are comparing a mult_g against a model.dual for the same solve, expect them to differ by a sign — that is by design, not a bug.

Before v0.9.0 the .sol writer emitted the Lagrange multiplier without converting it, so every AMPL/Pyomo dual came back negated relative to every other solver (#271). Objectives and primal solutions were never affected.

Exit codes

• 0 — Solve_Succeeded, Solved_To_Acceptable_Level, or Feasible_Point_Found. These are exactly the statuses whose solve_result_num falls in the 0–99 solved band, and the two channels are held to that by a test: the exit code and the .sol a run writes can never disagree about whether the solve succeeded.
• non-zero — any other ApplicationReturnStatus.

In AMPL solver mode (-AMPL) the exit code instead follows the AMPL contract: 0 for any solve that ran and produced a .sol file — limit-reached, infeasible, even a failed solve — since the termination is carried by the file’s solve_result_num. Genuine startup failures (unreadable .nl, bad option) still exit non-zero.

Diagnostics & introspection

pounce --about                                   # version, build info, features, backends
pounce problem.nl --dump kkt:5-10 --dump iterate # dump per-iteration diagnostics
pounce problem.nl --dump kkt --dump-dir /tmp/d   # override the dump root

• --about — print version, build info, enabled features, and linear-solver backends, then exit.
• --dump <cat>[:<spec>] — write the diagnostic category to per-iteration files (JSONL). Wired categories are kkt and iterate; an optional :<spec> selects iterations (e.g. kkt:5, kkt:2-10, iterate:all).
• --dump-dir <path> — override the dump root (default ./pounce-dump-<timestamp>).
• --dump-format <fmt> — dump format (default jsonl).

Help

pounce --help
pounce --version          # also -v, -V

Solver Options

POUNCE accepts options the same way upstream Ipopt does. Option names and semantics follow Ipopt’s, so an existing Ipopt options file or KEY=VALUE invocation works unchanged.

Setting options

On the command line — append KEY=VALUE pairs after the input:

pounce problem.nl tol=1e-10 max_iter=500 print_level=8

From an options file — upstream ipopt.opt format (one name value pair per line, # comments):

pounce problem.nl --options-file tuned.opt
pounce problem.nl option_file_name=tuned.opt   # the Ipopt spelling; same thing

From an options file nobody named. With neither of the above, POUNCE looks in the working directory for pounce.opt, then ipopt.opt, and reads the first one it finds — the way ipopt picks up ./ipopt.opt:

printf 'max_iter 5\n' > ipopt.opt
pounce problem.nl        # → Using option file "ipopt.opt".

The run says which file configured it (on the same sb gate as the banner). If both default names are present, pounce.opt wins and the other is named in a warning rather than left to look applied. --no-options-file skips the lookup entirely — the escape hatch for a directory holding an options file written for some other run.

Command-line KEY=VALUE pairs — and $pounce_options, AMPL’s <solver>_options channel — override values loaded from the options file, never the reverse.

An options file that was named but cannot be read is an error, not a shrug:

$ pounce problem.nl option_file_name=typo.opt
pounce: failed to load options file: options file "typo.opt" does not exist. …

Upstream opens a named file with a bare ifstream and reads nothing if that fails, so a typo there runs at stock defaults without a word. That silence is the one thing this path exists to remove (#518): a run configured entirely through an options file that quietly ran at defaults still reported success, which invalidates any benchmark set up that way.

One upstream quirk is inherited: option_file_name set inside an options file chains nowhere, because the file has already been chosen by the time it is read. POUNCE warns about that rather than ignoring it.

Commonly used options

For the full upstream option catalogue, see the Ipopt options reference; POUNCE reuses those names.

For scaling-specific options (nlp_scaling_method, target-gradient overrides, linear_system_scaling), see the Scaling reference page. For nonlinear bound tightening (presolve_fbbt, fbbt_tol, fbbt_max_iter, fbbt_max_constraints), see the FBBT reference page.

Options POUNCE does not implement

POUNCE’s option registry is a faithful port of Ipopt’s: every name Ipopt registers is registered here, so an ipopt.opt written for Ipopt parses unchanged. Registering an option is not the same as implementing it, though — and for a long time setting an unimplemented one did nothing at all, silently.

Options naming a feature POUNCE does not have now fail the solve, naming the option, the feature, and what to use instead:

$ pounce model.nl dependency_detector=mumps
pounce: `dependency_detector` configures linear-dependency detection on the
equality constraints, which pounce does not implement. It is registered so an
ipopt.opt written for Ipopt still parses, but setting it used to do nothing at
all — silently — so it is refused instead. Instead: pounce's presolve removes
structurally redundant rows; see `presolve`. Remove it to run.

The features in question: the Chen-Goldfarb (CG-penalty) / inexact-Newton line search, derivative approximation by finite differences, linear-dependency detection, the per-iteration NaN/Inf derivative check, multiplier recalculation by least squares, least-square initialization of all duals (least_square_init_duals), oracle-driven μ on the switch into fixed mode (fixed_mu_oracle; POUNCE implements that option’s default, average_compl), a selectable constraint-violation norm, magic steps, bound replacement, the L-BFGS augmented-system variants, skipping the finalize callback, the dynamic HSL loader, and suppress_all_output / debug_print_level.

Two line-search knobs joined that list with #551. theta_min is the CG-penalty acceptor’s threshold, not the filter’s — the filter line search derives its own theta_min from theta_min_fact * max(1, θ₀), as upstream does, and never takes it directly, so set theta_min_fact if that is what you meant. alpha_for_y_tol configures only the primal-and-full / dual-and-full multiplier-step rules, which POUNCE does not have; alpha_for_y supports primal (the default), bound-mult, min, max and full.

Four corrector knobs joined it too. corrector_type and its safeguards skip_corr_if_neg_curv, skip_corr_in_monotone_mode and corrector_compl_avrg_red_fact select and gate the corrector step Ipopt tries inside the line search (FilterLSAcceptor::TryCorrector). POUNCE’s line search takes no corrector trial at all. The predictor-corrector POUNCE does implement is Mehrotra’s, applied to the search-direction right-hand side: mehrotra_algorithm=yes is read and honoured, and it also selects mu_strategy=adaptive and mu_oracle=probing.

And three sub-capabilities of features that do run, where the refusal message has to say so or it reads as “restoration is missing”: expect_infeasible_problem_ctol / _ytol steer IpBacktrackingLineSearch’s count_successive_shortened_steps_ machinery, which POUNCE does not have — the restoration phase itself runs, and expect_infeasible_problem and required_infeasibility_reduction are read; resto_failure_feasibility_threshold asks for a threshold below which a stopped restoration is reclassified as a failure, and POUNCE has no such reclassification (max_resto_iter, below, is what bounds a restoration); and limited_memory_special_for_resto=yes asks for the special quasi-Newton update Ipopt dropped in Nov 2010 — L-BFGS runs in the restoration sub-solve with the regular update, which is what upstream’s own default no asks for.

The same rule applies one level down, to a single value of an option that otherwise works. bound_mult_init_method is read and honoured, but only its default constant is implemented; mu-based parses (an ipopt.opt written for Ipopt still loads) and is then refused, because serving it as constant would run a different initialization than the one you asked for under the name you asked for.

And it applies per entry point, for an option whose feature exists but sits somewhere this caller cannot reach. option_file_name is refused by a library caller that reads no options file, and the convex qp_* knobs (see LP/QP Routing) are refused by one that cannot route a model to the convex engines. Both are registered centrally, so they parse everywhere and the message can say which surface honours them — rather than the caller getting Unknown option for a name that exists.

option_file_name was on that list until #518 implemented it; refusing an option is the cheap half of “implement it or fail loudly”, and an entry leaves this table by getting the other half.

Two deliberate exceptions:

• Setting an option to its registered default is allowed. A generated ipopt.opt spells out defaults, and dependency_detector=none asks for nothing. Only a value that differs from the default is a request POUNCE cannot honour.

• Caching hints are checked, not trusted. grad_f_constant, hessian_constant, jac_c_constant and jac_d_constant tell the solver a derivative does not change between iterations. Ipopt takes such a hint on faith and silently returns a wrong answer if it is false. POUNCE asks the model first, and there are three cases:

  • POUNCE proves the derivative constant — from an .nl model’s own algebra — and reuses it across iterations whether or not you set the option. Setting it is harmless and unnecessary.
  • POUNCE proves the derivative is not constant and you set the option anyway: the option is ignored, with a warning. A QCQP’s ∇²L = σQ₀ + Σᵢλᵢ Qᵢ genuinely varies with the multipliers, so hessian_constant=yes there is not a hint but a false statement, and honouring it would trade a correct answer for a fast wrong one.
  • POUNCE cannot tell — every callback front end (the C interface, the Python Problem callbacks, both GAMS links) hands POUNCE numbers rather than algebra — and your assertion is honoured on trust, exactly as upstream. “Unproved” is not “disproved”; overriding you here would be its own silent wrong answer.

  POUNCE_DBG_CONSTDERIV=1 prints which of the three fired for each of the four options.

Options whose feature runs and whose value simply is not read yet are not in this category; they still solve, with the default in effect. Wiring those is tracked on #191 and #483.

Derivative checker

Wrong analytic derivatives are the most common cause of an NLP that stalls, cycles, or converges to something that is not a solution — and they are invisible from the iteration log. derivative_test compares what your TNLP returns against finite differences at the (bound-projected) starting point, before the solve:

pounce problem.nl derivative_test=first-order

first-order checks eval_grad_f and eval_jac_g; second-order adds eval_h; only-second-order checks the Hessian alone. The Hessian is checked one multiplier block at a time — obj_factor = 1, λ = 0 against differences of eval_grad_f, then obj_factor = 0, λ = eⱼ against differences of row j of eval_jac_g.

Entries that look wrong are marked *:

Derivative checker: first derivatives at the starting point (perturbation 1.0e-8, tolerance 1.0e-4).
* grad_f[    1]       =    3.5000000000000000e0    ~    3.0000000119209290e0  [  1.667e-1]
1 suspicious derivative(s) and 0 missing sparsity entrie(s) out of 6 checked (8 evaluations).

Two checks beyond upstream Ipopt’s, because both catch a class of bug no value-by-value comparison can:

• A Jacobian or Hessian entry whose finite difference is nonzero but which the sparsity structure omits (! in the report). A missing structural entry is not a wrong number — it is a derivative the solver can never see.
• The perturbation is taken downward when stepping up would leave a variable’s box, so a model using sqrt, log, or 1/x is not evaluated outside its own domain by the checker.

The test is advisory: it reports and the solve continues. It is written to stderr, so it survives print_level=0 and never mixes into --json-output’s stdout. It is slow — the second-order test costs roughly (m+1)·n evaluations — so leave it off for production runs.

check_derivatives_for_naninf is a separate upstream option, for a per-iteration NaN/Inf guard, and is not implemented.

Choosing a linear solver

POUNCE implements two KKT backends:

• feral — pure-Rust sparse symmetric indefinite solver. The effective default; no Fortran toolchain, no HSL licence.
• ma57 — HSL MA57, available only in a cargo build --features ma57 build.

The option’s registered value list is a faithful port of upstream Ipopt’s (ma27, ma77, ma86, ma97, mumps, pardiso, pardisomkl, spral, wsmp, custom), so an ipopt.opt written for Ipopt parses here unchanged. Selecting one of those fails the solve with a message naming it. They used to fall through to FERAL silently, which meant linear_solver=ma97 “worked” and a benchmark comparing backends compared FERAL with itself.

The registered default is feral, which diverges from upstream’s ma57 on purpose: a default has to name a solver the binary actually contains. Under the upstream default a pure-Rust build advertised MA57 to every print_user_options dump while running FERAL, and an HSL-enabled build used MA57 without being asked. If you build --features ma57 and want it, select it explicitly — that is the one behavioural change here.

Not a failure: explicit ma57 on a build without the feature falls back to FERAL and says so in the banner (FERAL (ma57 requested but not compiled)). That substitution is reported rather than hidden, and failing a portable ipopt.opt over a build flag would cost more than it buys.

The per-backend tuning options (ma97_scaling, mumps_pivtolmax, pardiso_*, wsmp_*, spral_*, …) remain registered for the same ipopt.opt-compatibility reason. They are unreachable now that their backend cannot be selected, so setting one alongside options POUNCE does read warns and solves — it does not fail the run:

$ pounce model.nl ma97_order=metis tol=1e-8
pounce: warning: `ma97_order` configures the HSL MA97 sparse symmetric linear
solver, which pounce does not implement, so it is ignored — as is every other
`ma97_*` option. pounce factors the KKT system with `feral` (pure Rust, the
default) or MA57 (`linear_solver=ma57`, in a `--features ma57` build); no
setting written for another backend transfers to either. The name is registered
so an `ipopt.opt` written for Ipopt still parses unchanged — which is why this
is a warning and not an error: the solve runs, and its result is unaffected.

A warning rather than a refusal, unlike the options above, because a portable ipopt.opt routinely carries settings for several backends at once so that one file runs everywhere. Refusing would fail that file over knobs the run never touches — for a user who is not using MA97 and never asked POUNCE to. One line is printed per backend family, listing the options it saw, and only for a value that differs from the registered default: a file that spells out defaults asks for nothing and gets nothing said about it. pardisolib warns with the pardiso_* family; hsllib is still refused, because POUNCE has an HSL backend (MA57) and the refusal points you at --features ma57 rather than leaving you to believe a library was loaded.

…unless they are all you set

That reasoning assumes the file has other business here. If the backend knobs are everything the run sets, nothing in the file survives, and warning-then-solving would answer “tune the linear solver” by tuning nothing and reporting success. That case is refused:

$ pounce model.nl ma97_order=metis
pounce: error: every option this run sets configures a linear-solver backend
pounce does not implement, so there is nothing left for it to act on. […] Set
`linear_solver=feral` (or `ma57`) if the defaults are what you want.

The rule in full:

Two details worth knowing. The second row counts an option’s presence, not whether you changed it — writing tol at its default is still a statement about this solve, and it is enough to put you back on the warning path. And option_file_name does not count as content: it says where the options came from, not what to solve, so pointing at a backend-only ipopt.opt is refused exactly as passing the same knobs on the command line would be.

A file that selects the backend it tunes never reaches this: linear_solver=ma97 is refused on its own, as above.

MA57 backend knobs

Only relevant in a --features ma57 build running linear_solver=ma57; see Installation. All of them are the upstream Ipopt names with the upstream meanings, so an ipopt.opt written for Ipopt-MA57 transfers unchanged.

Each can be scoped to the restoration sub-solve with a resto. prefix, e.g. resto.ma57_pivtol=0.5, which leaves the main solve’s value alone.

Up to and including 0.10.0 none of the nine reached the backend. They were registered, documented, parsed and validated — and then discarded: Ma57Options::from_options_list was a correct reader with no callers, because every construction site went through Ma57SolverInterface::new(), which hard-codes the defaults. Set to any value, the factorization behaved as though you had set nothing (#825). They are wired from 0.11.0 on, which means a build that was tuning MA57 through these knobs will now actually be tuned by them — expect the trajectory to move. ma57_pivtolmax also now refuses a value below ma57_pivtol rather than accepting an empty escalation range.

MA57 batched back-substitution (ma57_batched_backsolve)

Only relevant in a --features ma57 build running linear_solver=ma57.

Under the limited-memory quasi-Newton Hessian the solver builds a Sherman-Morrison-Woodbury correction from a handful of extra right-hand sides. Those could go to the linear solver in one call instead of one call per column, which saves traversals of the factor — but only if the backend’s multi-RHS answer is bit-identical to solving the columns one at a time, because these columns feed an iterate whose trajectory must not move. So the solver asks first, and a backend that blocks its triangular substitution across columns has to say no.

MA57 blocks. Its batched answer is tolerance-equal but not bit-equal, so it says no by default, and this option is how you overrule it.

Turning this on is a trajectory change, not a free speed-up. Measured on a 118276-row KKT system, the same binary with and without the batch diverges in the last printed digit of the objective at iteration 20 and finishes at a different iteration count. On a nonconvex problem a perturbation that size can select a different local optimum — issue #729 is MA57 taking pooling_rt2stp to an objective 25% worse while still reporting Optimal Solution Found.

What it buys, per iteration on that model:

against a 3–5% replicate spread, so the net is small and the headline row is not the whole story. Against Ipopt/MA57 on an equal-iteration basis it is more interesting: the back-solve row goes from 51% worse to 19% better.

Do not compare wall-clock across this option. The two settings walk different trajectories, so the difference is dominated by how many iterations each run happened to take rather than by work removed. In the measurement above the runs finished in 201 / 206 / 173 / 160 iterations, and the fastest arm was luck, not throughput.

Full write-up, including why the option has no width ceiling the way the FERAL backend’s equivalent does: dev-notes/ma57-batched-backsolve.md.

Withdrawing the constraint perturbation (perturb_delta_c_max_rungs)

POUNCE extension; not an upstream Ipopt option. Default 3; 0 restores the pre-#592 escalation exactly.

When the KKT factorization does not deliver the requested inertia, the solver climbs a ladder of perturbations: delta_w on the Hessian block, and delta_c on the constraint block. delta_c is the remedy for a rank-deficient constraint Jacobian, and it is reached for when the factorization reports Singular.

Since #540 a factorization also reports Singular when its inertia is unmeasurable — the count disagrees and the smallest pivot sits at the noise floor. That is evidence about the measurement, not about the Jacobian’s rank. When the Jacobian in fact has full rank delta_c cannot help, and because it stays switched on for the rest of that augmented system, the delta_w ladder then has to climb against a matrix delta_c has made harder to hit the requested inertia on. On the #592 model that cost five rungs, ending at delta_w = 1e2 where Ipopt accepted the step at 1e-4; the over-damped step froze the objective for eight iterations and the solver exited at a point a restart improved by 0.08%.

Rather than predict which kind of Singular a report was — the counts are the very thing #540 established are noise — the ladder answers it empirically. After this many rungs with delta_c on and still no acceptable inertia, delta_c is withdrawn, the delta_w ladder restarts, and delta_c is latched off for the remainder of that augmented system; the next iterate starts clean. Lower values withdraw sooner.

Where delta_c is the right remedy this never fires: on eigena2 and eigenb2 it is followed by at most one rung. See crates/pounce-common/src/pd_perturbation.rs (maybe_withdraw_delta_c) and dev-notes/issue-592-restart-non-idempotence.md.

Inertia-free curvature test (neg_curv_test_tol)

By default every KKT factorization is checked for the right inertia — as many negative eigenvalues as there are constraints — and the primal regularization δ_x is escalated until it has it. The inertia-free alternative of Zavala & Chiang (2014) factors without that check and instead asks whether the direction the system produced actually curves upward:

dxᵀ W dx + dxᵀ Σ_x dx + dsᵀ Σ_s ds [+ δ_x‖dx‖² + δ_s‖ds‖²]
    ≥ neg_curv_test_tol · (‖dx‖² + ‖ds‖²)

This is a heuristic, and turning it on is not free — it can change the answer, not just the path to it. Measured over POUNCE’s fixture corpus at the recommended 1e-11 (scripts/sweep-fixtures.sh, both legs), 11 of 59 models move:

The last four rows are the reason to read this before switching it on. deb7 and eigenb2-under-L-BFGS stop converging at all; autocorr_bern55-06 and pooling_rt2stp still report success but land on a worse objective — a tolerance-legal wrong answer, which is the failure mode that is invisible to a suite asserting status and objective-to-a-tolerance. Accepting a factorization whose inertia is wrong is exactly the kind of change that produces it.

It is off by default (neg_curv_test_tol=0 keeps the inertia check), and nothing above happens to a solve that leaves it alone. If you turn it on, measure your own model.

Escaping a stationary point that is not a minimum (neg_curv_escapes)

The convergence test is a first-order test. On a nonconvex model that is strictly weaker than “local minimum”: at a point where the reduced Hessian on null(A) is negative definite every KKT residual is zero, so the test has nothing to object to, and the point reported as Solve_Succeeded can be a constrained maximum.

The CLI fixture nonconvex_qp.nl is that case in three lines:

min x₀·x₁   s.t.   x₀ + x₁ = 2,   0 ≤ x ≤ 4

On the feasible segment the objective is f(x₀) = x₀(2 − x₀), which is concave — maximized at (1, 1) with f = 1, minimized at the endpoints (0, 2) and (2, 0) with f = 0. From the bound-pushed start (0.01, 0.01) the first Newton step lands exactly on (1, 1), and every iteration after it takes a step of size 1e-14.

Inertia correction does not prevent this and never could. It engages — the iteration log shows lg(rg) from the second iteration on — but δ_x I is symmetric, the model and the iterate are symmetric under x₀ ↔ x₁, and a symmetric correction applied to a zero gradient gives a zero step however indefinite the reduced Hessian is. Regularization makes the step well-posed; nothing else in the algorithm asks whether the point it converged to is a minimum.

With this on, a point about to be certified is first tested for second-order necessity: one extra factorization of the augmented system with the inertia check on and no perturbation, whose correct inertia is exactly the statement that W + Σ is positive definite on null(A). A point that passes costs that one factorization and nothing else. A point that fails gets δ_x escalated until the inertia is right, and a few inverse-iteration back-solves against that factor recover the most-negative-curvature direction — which is then measured, not trusted. The solve steps along it (capped by the fraction-to-the-boundary rule, backtracked against the second-order decrease model, refused outright if it raises the constraint violation past constr_viol_tol) and continues.

It cannot return a worse answer than leaving it off would have. The stationary point is snapshotted before the step and is restored and reported unless the continuation comes back with a certificate of its own at a better point — the same floor-and-deadline accounting as resto_decline_deferrals, and the continuation is cut after 30 iterations either way. Raising the option above 1 does not weaken that (gh #805): the floor holds the best certificate the escapes have left, not the most recent one, so every bet is placed against the same baseline the first one was — the point a neg_curv_escapes = 0 build reports. Each escape does buy its continuation its own 30 iterations, so the cost scales with the option and the guarantee does not. On nonconvex_qp.nl — and on nonconvex_qp_ineq.nl, the same model with its row relaxed to x₀ + x₁ ≥ 2 — it turns Solve_Succeeded at obj = 1 into Solve_Succeeded at obj = 0; across the rest of the fixture corpus (scripts/sweep-fixtures.sh, both legs, 152 fixture-legs) it moves nothing.

Two limits are worth knowing:

• It is still a local method. An escape finds a point that is second-order suspect and leaves it; it does not certify global optimality, and a stationary point whose reduced Hessian is positive definite is never touched.
• Under hessian_approximation=limited-memory it does nothing. The curvature it reads is B, and BFGS maintains B positive definite by construction, so the inertia test passes at δ_x = 0 and the escape declines. The L-BFGS leg of the fixture sweep still reports obj = 1 on both nonconvex-QP fixtures.

Bound relaxation and honor_original_bounds

Before the solve, the NLP arm widens every variable and constraint bound by bound_relax_factor (default 1e-8, capped by constr_viol_tol), exactly as upstream Ipopt does — a feasible-iterate log-barrier needs x strictly inside its bounds, and this keeps the iterates there without the user’s bounds becoming numerically degenerate. The consequence is that a solution pinned to a bound is reported just past it:

min (x − 3)²  s.t.  0 ≤ x ≤ 1     →     x = 1.00000000937

honor_original_bounds=yes projects the reported point back into the bounds you declared, so that solve returns exactly x = 1. Reach for it whenever the value flows somewhere that cares about the domain — a sqrt(1 − x), a domain assertion, or a Pyomo Var the value is loaded back into.

The default is no, matching upstream. As upstream also documents, the constraint-violation and complementarity figures in the end-of-run summary are for the non-projected point; only the reported x (and the objective and constraint values evaluated at it) move.

Note what honor_original_bounds does not do: it moves the reported point, not the solve. The iterate still stopped where the relaxed bounds put it, and any question of the form “is this constraint active” is still being asked about a point ~1e-8 shy of the bound. Projection cannot recover that, because the projection has no way to tell “pinned to the bound” from “genuinely 1e-8 inside it”. Crossover answers that question instead: it re-solves against the bounds you declared, so the returned point sits on them and the active set is established rather than inferred.

The convex arm does not widen by default

The LP / convex-QP / conic arms solve the model exactly as declared. Widening moves the optimum by δ times the bound’s multiplier, and nothing bounds that product: on LISWET1 — every one of 10 000 monotonicity rows active, multipliers summing to 1.6e9 — a 1e-8 widening buys 9.0 of objective, a 33 % error against the published optimum. Scored on the Maros–Mészáros optima (DOC 97/6) the convex arm is 138/138 correct without it and 130/138 with it, so it is off by default there.

Set bound_relax_factor explicitly and the convex arm applies it and reproduces the NLP arm’s model exactly, so Ipopt parity remains available on request.

The two arms therefore disagree on constraint-degenerate models, by design. final_declared_constr_viol reports how far outside the model as declared a returned point sits, on either arm, so the difference is readable rather than silent; the console prints it as one extra line whenever it differs materially from Constraint violation.

A convex solve that declines to certify is handed to the NLP arm (gh #535), and that re-solve also runs on the declared model unless you asked for a widening by name.

Large constraint values and primal_noise_floor_kappa

On a model whose constraint values run to ~1e7 and beyond, a converged solve could exit Search_Direction_Becomes_Too_Small while holding the correct optimum. The cause is arithmetic, not the model.

The KKT error the convergence test compares against tol is

max( ‖∇L‖∞ / s_d ,  max(‖c‖∞, ‖d − s‖∞) ,  ‖compl‖∞ / s_c )

The dual and complementarity terms are normalised; the primal one — like upstream Ipopt’s — is a bare absolute residual. But c_i = g_i(x) − b_i and d_i − s_i are each a difference of quantities the row’s own size, so they are quantised in units of eps · |b_i|. At |b| ~ 1e8 the smallest nonzero value the primal term can take is one ulp, 1.5e-8 — already larger than the default tol = 1e-8. Asking for nlp_err <= tol there is asking the residual to land on a bitwise-exact 0 rather than on one ulp, which is arithmetic luck rather than a property of the iterate.

POUNCE therefore judges the primal term in the strict test against each row’s own floating-point resolution: a row’s residual counts only where it exceeds max(placement floor, kappa · eps · |row magnitude|), with kappa = primal_noise_floor_kappa (default 64). Verdicts are flat across kappa from 8 to 1024 on the measured set.

Three things bound what this can do:

• Only the strict test reads it. constr_viol is still checked against constr_viol_tol (default 1e-4) on the full, unfloored residual, so nothing the floor forgives can exceed the feasibility tolerance you set — however large your data grows.
• The acceptable-level band keeps the raw error. It sits two decades above tol, clear of any realistic quantum.
• It cannot rescue an infeasible model. On a model with no feasible point the filter and restoration phase reach a verdict on their own criteria; the floor only ever participates at a point the rest of the algorithm already believes is converged.

Set primal_noise_floor_kappa = 0 to switch the floor off and restore upstream Ipopt’s bare-absolute primal term exactly.

When the floor changes the reported picture, the end-of-run summary says so — a large-|b| solve prints the tested value under the raw one:

Overall NLP error.......:   2.3841857910156250e-07    2.3841857910156250e-07
  ...above the per-row floating-point noise floor:   0.0000000000000000e+00

Solves where the two agree — every model whose data is O(1) — print the usual block unchanged.

One case this does not paper over: tightening constr_viol_tol below a row’s own ulp (say 1e-8 on data at 1e8) still will not certify. That is the tolerance gate doing what you asked — the residual you requested is not representable at that scale.

Solved_To_Acceptable_Level and acceptable_progress_kappa

Solved_To_Acceptable_Level is the fallback verdict for a solve that cannot reach tol: after acceptable_iter (default 15) consecutive iterates with an NLP error under acceptable_tol (default 1e-6), the solver stops and hands back the point it has. That criterion is a count of iterates inside a band, and on its own it asks only is the error small — never has anything stopped moving.

Those come apart. An interior-point iterate can be near-stationary for the current barrier subproblem — a much weaker statement than near-KKT for the NLP — for fifteen iterations running while the solve is still descending. Two measured cases: the kissing model stopped with objective 1.00000108 where continuing reaches 0.84544259 and a strict certificate, 18% lower; NARX_CFy stopped with both residuals near 1e-7 where sixty more iterations collapse them by five orders.

POUNCE therefore also requires the streak to have flattened. Across the acceptable_iter iterates that made it up:

• the spread (max − min) of the NLP error must be within acceptable_progress_kappa · acceptable_tol; and
• the spread of the objective within the same fraction of acceptable_tol · max(1, |f|).

acceptable_progress_kappa defaults to 0.1, so at default tolerances both quantities must have stayed inside a tenth of the acceptable band over the whole streak.

It is a spread, not a trend, and either signal alone is enough to keep solving. Both choices are deliberate: kissing’s error was an order of magnitude worse at the iterate it stopped on than at one it had already reached inside the same streak — it was wandering across the band, not converging inside it — while its objective was flat to all eight printed figures over the same iterates.

Three things bound what this can do:

• It cannot lose a verdict. The refused termination is recorded, and a run that fails to do better ends at exactly that iterate under exactly that status. A misfire costs iterations, never the answer — you will not see Maximum_Iterations_Exceeded where the count alone would have said Solved_To_Acceptable_Level.
• It never looks at a solve that converges. A solve that reaches tol never completes an acceptable-level streak, so nothing here runs.
• A genuine stall flattens. When the iterate, the objective and the error are all pinned — the case the acceptable-level exit exists for — the window is flat and termination happens as before.

Set acceptable_progress_kappa = 0 to switch the progress test off and restore upstream Ipopt’s bare consecutive-count criterion. Widening acceptable_tol widens the flat bar with it, so asking for a looser band still gets you the early exit.

A settled point with a runaway multiplier (dual_divergence_retry)

Some models have a solution at which a constraint’s gradient vanishes. The row is satisfied, the primal iterate is exact — but the multiplier that would certify it is arbitrary rather than nonexistent, and the barrier drives it off to infinity.

The standard case is a complementarity constraint lowered as a product, G(x)·H(x) = 0, at a point where the pair is biactive: G = 0 and H = 0 together. The product’s gradient is H∇G + G∇H, and both terms vanish there. MPCC lowerings (ncp_eq, prod_eq) reach such points routinely; MacMPEC’s qpec_small does at its solution (1, 1, 0).

This used to ship a wrong verdict, because the convergence gate reads an NLP error normalised by s_d, and s_d grows with the mean multiplier magnitude — so a runaway multiplier divides itself out of the number the gate tests. One line of the summary block, on qpec_small at bound_relax_factor=0:

                                   (scaled)                 (unscaled)
Overall NLP error.......:   8.2335532426389998e-11    7.8965510781517834e+04

Fifteen orders apart. The gate read the left column and reported Solved_To_Acceptable_Level.

POUNCE now watches for the signature directly. At one and the same iterate:

• the primal infeasibility is at zero (≤ 1e-8);
• the step has settled — maxᵢ |dᵢ| / (1 + |xᵢ|), over the x and s blocks, is at or below dual_divergence_retry_step_tol (default 1e-5); and
• the unscaled Lagrangian-gradient norm is at or above dual_divergence_retry_du_floor (default 1e2).

The middle conjunct is the one doing the work, and it is what separates this from an iterate that is simply diverging: a diverging solve has a large step, not a zero one. It is also the only barrier protecting models with no sign-feasible multiplier, where the remedy below reaches a plausible-looking answer below the true optimum. Measured: qpec_small settles to 4.3e-8, MacMPEC ralph1 — which has no multiplier at all — bottoms out at 7.2e-3. The 1e-5 default sits between them, and a fixture pinning that gap ships in crates/pounce-algorithm/tests/issue_884_biactive_dual_divergence.rs.

Only a model with at least one constraint row is eligible: on an unconstrained model ∇L ≡ ∇f, so the third conjunct would be a second, much looser copy of dual_inf_tol.

When the signature is seen and the solve ends Solved_To_Acceptable_Level or Restoration_Failed, POUNCE re-solves once from cold with perturb_always_cd=yes — regularising the constraint block from the first factorisation rather than waiting for an inertia failure that a vanishing gradient never produces. The second answer is returned only if all of:

1. the retry ends Solve_Succeeded;
2. its unscaled KKT error and constraint violation are both within acceptable_tol;
3. its unscaled KKT error is strictly better than the first attempt’s; and
4. its answer is admissible next to the first attempt’s — it may not return a strictly worse objective, and an objective improvement may not be bought with primal slack (a better objective at a larger constraint violation is refused). Both comparisons use acceptable_tol scaled by max(1, |first objective|), and both stand down when the first attempt is not itself feasible within that tolerance, since there is then no admissible answer to protect.

Those two statuses are the whole scope, and they are the two a vanishing gradient row produces directly: Solved_To_Acceptable_Level is gh#884 verbatim, and Restoration_Failed is the same defect one step earlier. Error_In_Step_Computation and Maximum_Iterations_Exceeded are deliberately not on the list even though the detector can legitimately fire before them — they are generic exhaustion exits that any hard model can reach for unrelated reasons, and retrying there buys nothing while costing a full second budget. See the cost note below for the measurement that set this.

Otherwise the first attempt’s status, point, statistics and final trace row are all put back — the retry costs iterations, never the answer. Condition 2 is what stops the gate reproducing the bug one attempt later: the defect was a status its own unscaled residual contradicted, so a promotion rule reading the status alone would launder it again.

Condition 4 is a separate barrier, not a restatement of the others. Conditions 1–3 rank the two attempts on their certificates, and a certificate cannot say which of two feasible points you should receive: any other KKT point satisfies the KKT conditions in the model’s own units just as well. Without condition 4, measured over 400 random QPECs under the exact-product lowering, 42 of 68 promotions returned a different local solution and three returned a strictly worse feasible point — worst case -13.0057 given up for -1.2072. On MacMPEC’s scholtes4, whose optimum is exactly 0, the retry returned -6.61e-05 — unreachable by any feasible point — by moving the complementarity row from 2.07e-25 to 1.09e-09, and reported Optimal Solution Found.

When condition 4 is what refuses a retry, the console says so explicitly — “declined on the ANSWER, not the certificate” — with both attempts’ objectives and constraint violations, because a converged retry with a clean certificate being refused otherwise reads as a contradiction.

On qpec_small the retry takes the unscaled KKT error from 7.8966e+04 to 9.9636e-08 — nine orders — at the cost of a primal residual that goes from 1.1e-16 to 5.5e-12, and a point 3.7e-06 further from (1, 1, 0). That trade is the point: a marginally looser answer that comes with a certificate a reader can check.

Cost. The detector runs once per iteration and reads quantities the convergence check already computes. The retry itself is the outermost wrapper, so it never runs where an inner one already won — at default options qpec_small is rescued by the μ-strategy fallback (see “Barrier-parameter (μ) strategy” below) and the dual-divergence retry spends nothing. Worst case is one extra solve, under your own max_iter, on a run that was already reporting a non-success verdict.

Among acceptable-level exits in the 80-fixture regression corpus, on both sweep legs, nothing else reaches the floor: the closest non-MPCC approach is eigena2 under L-BFGS at an unscaled dual of 37, and its step is 7.9e-9 — settled, but two orders under the 1e2 floor. One fixture does reach the floor at another status, and it is why the scope names statuses: deb7 under L-BFGS settles at iteration 346 to a step of 6.5e-6 with an unscaled dual of 9.2e+05 — above qpec_small’s on the dual conjunct, so no floor excludes it, and separable on the step conjunct only by tightening the default onto one fixture and spending the margin that holds ralph1 out. There the detector is right and the remedy is not: an earlier build that retried on Error_In_Step_Computation spent 715 → 3000 iterations to return the same status and the same objective.

Scoping by status is only as complete as the status is stable, and on that same fixture it is not: under limited_memory_ls_failure_restarts=1 (off by default) deb7 exits Restoration_Failed instead of Error_In_Step_Computation, so it is in scope.

So there is a second gate, and it reads the answer rather than the trajectory. The detector fires on an iterate; nothing in it says the solve ends there. A run can pass through a settled point with a diverged multiplier, work its way back down, and report something ordinary — and then there is nothing left for perturb_always_cd to repair.

What #884’s defect looks like in the answer is a point converged except that one multiplier ran away: the primal is exact, complementarity is met, and the whole residual is dual infeasibility. So the retry runs only when the reported answer’s unscaled constraint violation and unscaled complementarity are both at or below 1e-6 times its unscaled dual infeasibility:

Twelve orders. deb7’s complementarity is five percent of its own KKT error — that answer is not a converged point with a runaway multiplier, it is an unconverged point, and before this gate that run paid a full cold re-solve (6.1 s to 25.2 s) to decline an answer that was never going to be promoted (#887).

The test is a ratio within one answer on purpose. A floor on the reported residual would be a threshold on a scale-dependent quantity, and it does not even separate these cases — deb7’s 9.9e+01 sits one percent under the detector’s own 1e2. Comparing the answer against the runaway the detector saw does separate them, but it reads two numbers from a trajectory, and a trajectory is not stable across platforms. A ratio between two residuals of the same answer carries no units and cannot move that way.

One consequence is worth knowing if you are reading a report from a hard model: whether this gate opens is a property of the answer, not of the machine, but which answer a hard model reaches can differ between platforms. deb7 under that rung is the measured example — objective 99.677 on macOS against 99.651 on Linux, and on Linux the answer it reaches genuinely does carry the runaway, so the retry runs there and is supposed to.

Set dual_divergence_retry=no if you are running a hard model to a failure verdict and even one extra attempt is not worth the clock.

Two off switches, coarse and fine:

• dual_divergence_retry=no disables the retry outright. The detector still runs and the report still records what it saw.
• dual_divergence_retry_step_tol=0 holds the detector off, so nothing downstream of it can fire.

Both the signature and the retry’s outcome are reported. The summary block prints

Biactive dual divergence (gh#884)                    = detected

when the signature was seen and the solve did not end Solve_Succeeded (passing through such an iterate and recovering is routine on an MPCC, and a warning over a correct answer is noise), and the retry prints its own verdict line when it runs. In the JSON report the two are statistics.dual_divergence_signature and statistics.dual_divergence_retry_promoted.

Big models that start feasible and the theta_max ceiling

The filter has a hard ceiling. Any trial iterate whose constraint violation θ exceeds

theta_max = theta_max_fact · max(1, θ₀)

is rejected outright, before any of the filter’s usual tests run. It is a global-convergence safeguard: it keeps the line search from wandering arbitrarily far from feasibility.

The trouble is the 1. POUNCE’s θ is a 1-norm over constraint rows — ‖c‖₁ + ‖d − s‖₁, a sum of m residuals — so a ceiling of T really says “a mean per-row violation of T/m”, and that allowance shrinks as the model grows. And on a problem started at a feasible point, θ₀ = 0, the max collapses and the ceiling is the bare constant theta_max_fact however large the model is.

robot_a is the measured case: 52 013 constraint rows, a feasible start, so theta_max locked at 1e4 — a mean per-row allowance of 0.19 — while the route to the optimum passes through θ ≈ 9.4e7. Every step toward the solution was refused at the gate, and the solve ground to its iteration limit at objective 8.173304 instead of the true 1.0431952.

POUNCE’s answer is theta_max_adaptive_trigger (default 3), described below. theta_max_row_scale_kappa (default 0, off) is an earlier, static attempt at the same problem, kept because it is occasionally the more direct lever; it floors the reference at the row count instead:

theta_max = theta_max_fact · max(θ₀, theta_max_row_scale_kappa · rows, 1)

so the ceiling means a mean per-row violation of theta_max_fact regardless of m. Measured under defaults otherwise, against Ipopt 3.14 on the same machine:

Ipopt has the same defect and no correction for it; all three solve under theta_max_fact = 1e8 set by hand, which is the blunt version of the same move.

The adaptive rule: theta_max_adaptive_trigger

On by default. Rather than guessing from problem size whether a model needs headroom, POUNCE measures whether the ceiling is what is refusing the line search, and raises it only then.

A trial refused because θ_trial > theta_max takes a distinct early exit, before the filter and Armijo tests run at all. So the acceptor can count those refusals and compare them against the number of trials attempted. When every trial of a line search was refused at the gate, for theta_max_adaptive_trigger consecutive line searches, the ceiling is demonstrably the binding constraint — not the filter — and it is multiplied by theta_max_adaptive_factor (default 100), at most theta_max_adaptive_max_raises times per solve (default 4).

Three properties follow, and they are what the static floor could not offer:

• A converging model cannot trip it. Converging means trials are getting past the gate; a model accepting steps never accumulates the streak. brainpc1/3/5/7 — the family the static floor damaged at every kappa — are untouched by construction, not by a lucky constant.
• A blocked model trips it immediately. robot_a is refused at the gate from its first line search onward.
• The ceiling stays finite. Wächter–Biegler’s global-convergence argument (Thm. 2) needs theta_max finite, not fixed. A bounded number of bounded raises keeps it finite, so a solve cannot ratchet the safeguard away one line search at a time.

Requiring a streak rather than a single line search is deliberate: one Newton direction that overshoots into a huge θ can legitimately have all its trials refused, and backtracking is the right response to that. Only a model that cannot get past the gate repeatedly is one whose route needs the headroom.

Measured, defaults otherwise:

Across the whole Vanderbei corpus (733 problems) the rule changes four outcomes: britgas goes from its iteration limit to Optimal in 16 iterations, catenary solves to the same objective in 50 iterations instead of 56, and coshfun and brainpc0 fail either way — coshfun now reporting diverging iterates, which is what Ipopt 3.14 also does on it. Net Optimal count 702 → 703.

Note brainpc0 does trip the rule while brainpc1/3/5/7 do not, despite identical row counts. That is precisely the distinction a size-based floor cannot draw.

The restoration sub-IPM always runs with the rule disabled. Upstream already corrects the resto phase’s instance of this degeneracy by hard-coding resto.theta_max_fact = 1e8 (IpRestoMinC_1Nrm.cpp:91), so a rule that ratchets further would be compounding a correction already made.

Set theta_max_adaptive_trigger = 0 to restore upstream Ipopt’s fixed ceiling exactly.

When to reach for the static floor instead

Symptoms, all three together:

• a large number of constraint rows (thousands upward);
• a feasible or near-feasible starting point — the iteration log’s first inf_pr is 0 or very small;
• the solve stalls with inf_pr flat and alpha_pr tiny, and raising max_iter does not help.

The quick confirmation is to set theta_max_fact = 1e8 by hand. If that unsticks the model, theta_max_row_scale_kappa = 1 is the principled version of it — it scales the ceiling to the model rather than to a constant you picked.

Why the static floor is off by default

Because raising the ceiling unconditionally is not free. It relaxes a global-convergence safeguard, and a model that was not being blocked by it can wander instead. On the Vanderbei corpus, brainpc1/3/5/7 (m = 6900, θ₀ = 1e-2) all regress — brainpc1 from Optimal in 64 iterations to divergent, objective 3.7e3 against the correct 4.4e-04.

A scan over kappa showed why this cannot be tuned away. The damage is a step function, not a gradient:

brainpc3 and brainpc7 land on the identical worse answer at every nonzero kappa, even 0.01 — where the ceiling moves only from 1e4 to 6.9e5. The instant it rises at all, they break. robot_a meanwhile improves monotonically all the way to kappa = 1. There is no separating value.

That is a verdict on the design, not on the tuning: the real question is whether a model’s route to the optimum needs the extra headroom, and the row count does not answer it. A static floor cannot know — which is why the adaptive rule above, which asks the question directly, is the default and this one is not.

What still bounds the option when you do turn it on:

• It only ever raises the reference. A model whose own θ₀ already exceeds kappa · rows gets exactly upstream’s ceiling.
• It reduces to upstream on a single-row problem, where the floor is max(kappa · 1, 1) = 1 — upstream’s constant.
• theta_max is still finite, and still fixed for the whole solve after its first line search. This rescales the safeguard; it does not remove it.
• The restoration sub-IPM is untouched, at any kappa. Upstream already fixes its own instance of this by hard-coding resto.theta_max_fact = 1e8 (IpRestoMinC_1Nrm.cpp:91) — the resto NLP is also initialised feasible, so it hit the same degeneracy. Stacking the row floor on top would push that inner ceiling to 1e8 · m, i.e. remove it, so the sub-IPM always runs with kappa = 0.

Large gradients and dual_inf_scale_kappa

The dual side of the same story. dual_inf_tol (default 1.0) is a bare absolute bound on ‖∇L‖∞ — but the aggregate above normalises that quantity, dividing it by s_d, which grows with the mean magnitude of the multipliers. On a model whose gradients live at 1e10 the two gates are judging one number by standards ten orders apart.

Vanderbei’s orthrds2 is the reported case: s_d ≈ 1.6e10 with ‖∇L‖∞ = 89.7, so the aggregate’s dual term is 5.6e-09 — comfortably inside the default tol = 1e-8, i.e. stationary to nine digits relative to the size of the gradients involved — while the component gate refused it against 1.0. The solve exited Solved_To_Acceptable_Level holding the answer, and dual_inf_tol=1e3 alone turned it into Optimal Solution Found at the same objective.

The simplest statement of the defect: multiply an objective by a positive constant. Same feasible set, same solution, same active set, same Newton step — and every multiplier, s_d and ‖∇L‖∞ scale with it, so a large enough constant costs the certificate.

The strict test therefore judges the unscaled dual infeasibility against

max( dual_inf_tol ,  kappa · tol · dual_scale )

with kappa = dual_inf_scale_kappa (default 1) and dual_scale the magnitude of the largest single term ∇L is assembled from (∇f, Jᵀy, the bound multipliers). Since ∇L is the sum of those terms, ‖∇L‖∞ / dual_scale is the fraction of them that failed to cancel — a scale-invariant statement of stationarity, and the thing the absolute bound was standing in for.

What bounds it:

• It cannot forgive non-stationarity. A point where nothing cancelled has ‖∇L‖∞ ≈ dual_scale, a ratio of 1 against a bar of 1e-8. min −exp(x) s.t. x >= 0 running away to inf_du = 8.8e+47 is refused by eight orders, because its ∇f runs away by exactly the same factor.
• The aggregate still has to pass. nlp_err <= tol is tested on the same iterate; this only removes the second, inconsistent standard.
• It is inert on ordinary models. At the defaults the floor does not rise above dual_inf_tol until dual_scale exceeds dual_inf_tol / tol = 1e8, so every model with O(1) gradients keeps upstream’s comparison bit for bit.
• Only the strict gate reads it. acceptable_dual_inf_tol (1e10) is untouched.

Set dual_inf_scale_kappa = 0 to switch the floor off and restore upstream Ipopt’s bare-absolute bound. That is also the setting to reach for if you tighten dual_inf_tol and want that absolute standard honoured unconditionally — the floor is a floor, so it can override a tightened dual_inf_tol on a large-gradient model.

s_max — where s_d and s_c come from

The two normalising factors above are built from the multipliers themselves, capped by s_max (default 100, upstream’s):

s_d = max( s_max , (‖y_c‖₁+‖y_d‖₁+‖z‖₁+‖v‖₁) / (their total dimension) ) / s_max
s_c = max( s_max , (‖z‖₁+‖v‖₁) / (their dimension) ) / s_max

Both are exactly 1 while the multipliers average below the cap — which is every well-scaled problem, and why the option is invisible there — and grow as mean / s_max once the average passes it. Raising s_max therefore delays the normalisation (the KKT error stays closer to the raw residuals); lowering it applies the normalisation sooner and makes the scaled error smaller for the same iterate, so the solve certifies earlier. The scaled and unscaled numbers are both reported: --json-output carries final_kkt_error and final_unscaled_kkt_error, and their ratio is exactly what s_max controls.

Objective sense and obj_scaling_factor

obj_scaling_factor multiplies the objective the IPM minimizes, so a negative value maximizes — upstream’s documented spelling for a maximization problem stated as a minimization. Because it changes what is being optimized rather than just its conditioning, it is honored only by the general NLP interior-point path: a model that would otherwise route to the specialized convex solvers (LP / convex QP / SOCP, see LP/QP Routing) is re-routed under solver_selection=auto, and an explicit convex solver_selection is refused rather than silently answering with the minimizer.

A positive factor is a pure conditioning knob; the convex path reports natural units either way, so it keeps the fast path.

Starting-point conditioning

Three options displace the starting point before the barrier solve, and one turns the automatic retry that uses them on or off. All are off by default except the retry, which fires only after a solve has already failed. Full rationale and the measurements behind the defaults: Conditioning the starting point.

start_point_conditioner=adam reads three more. They are ignored unless it is set, and their defaults are KRONOS’s published stage-0 values (Ahmed & Hasan 2026, see Acknowledgments).

The warm-up is guarded: if it does not reduce the merit it hands back the original point unchanged. It is off by default because across 40 problems POUNCE already solves it cut the median iteration count to 0.83× while raising the total 1.62×, on the strength of two outliers (palmer1c 71 → 1023). A median win with a 14× tail is an option, not a default.

Barrier-parameter (μ) strategy

The barrier parameter μ controls the inner subproblem’s relaxation of complementarity. The two strategies are monotone (default — geometric schedule) and adaptive (quality-function oracle picks each μ from the current iterate’s complementarity). See μ-strategy for when to switch.

Quality-function oracle (adaptive-μ details)

These are only consumed when mu_strategy=adaptive and mu_oracle=quality-function. Defaults mirror upstream IpQualityFunctionMuOracle::RegisterOptions.

Adaptive-μ globalization

Tuning the safeguards that fall back to monotone-μ mode when the adaptive oracle stops making progress. Defaults mirror upstream IpAdaptiveMuUpdate::RegisterOptions.

Hessian approximation (hessian_approximation)

Which second-derivative information the algorithm works from. Four values; the first two are Ipopt’s, the last two are POUNCE extensions with no upstream counterpart.

exact and limited-memory are documented under L-BFGS initialization and throughout this page. The other two are below.

finite-difference

Recovers the exact Lagrangian Hessian by differencing the analytic Jacobian along a set of probe directions, rather than approximating it from step/gradient pairs. Where L-BFGS presents the linear solver a dense low-rank correction over a diagonal, this hands it the model’s real sparse Hessian, so a structured problem factors like one.

The cost is one Jacobian evaluation per probe group, per rebuild. Three options control how many groups there are and how often you pay for them:

fd_hessian_pattern. declared uses the TNLP’s declared Hessian structure — the structure call only, never the values, so it is available to any model that cannot evaluate second derivatives; every .nl declares one through AMPL’s AD. jacobian derives it from the Jacobian pattern alone, as the union over rows of supp(∇g) ⊗ supp(∇g). That is a strict superset of the true pattern, so it is safe — a superset costs extra probe groups, never a wrong answer — but it is not free: on benchmarks/large_scale laptime it is 146 267 nonzeros against the true 28 000. There is deliberately no mode that guesses a subset, which would silently drop curvature.

fd_hessian_coloring. A star colouring needs fewer groups than Curtis-Powell-Reid — on the Jacobian-derived laptime pattern, 42 where CPR needs 76 — and its recovery is algebraically exact, so it looks like the better default and is not. A forward difference is not an exact Hessian-vector product: it carries a third-derivative cross term into each row, and CPR’s distance-2 property forbids the two columns that term needs from sharing a group, while a star colouring does not. Measured on laptime, star over the Jacobian pattern takes 404 iterations to a wrong objective where CPR takes 38 to the right one; over the sparser declared pattern both take 30. Group size is not the cause — declared/star has the largest groups of the four and is fine. Use star only on a sparse declared pattern, and measure.

fd_hessian_reuse_tol. A rebuild costs one Jacobian evaluation per group, so skipping it when nothing has moved is the cheapest saving available. Both the primal iterate and the multipliers are tested, not just x: the Lagrangian Hessian is ∇²f + Σⱼ yⱼ ∇²cⱼ, so a cached Hessian is stale the moment y moves even if x has not — and the endgame of an interior-point solve is full of short steps with moving duals. 0, the default, rebuilds every iteration.

finite-difference will produce a Hessian for a model whose second derivatives do not exist, because it never asks for them. hessian_approximation=exact refuses such a model; this one does not. That is the point of it, and also the risk: the answer is only as good as the Jacobian is differentiable.

partitioned

Keeps one small dense quasi-Newton block per element function — the objective, and each constraint row, whose support is a row of the Jacobian — and assembles them into a genuine sparse Hessian. The linear solver then sees the model’s real block structure instead of the diagonal the limited-memory low-rank path presents. Intended for structured problems where second derivatives are unavailable but the Jacobian sparsity is declared; direct-collocation trajectory optimization above all.

partitioned_elements. per-constraint gives each block a multiplier-independent target and assumes nothing about variable ordering, at the cost of as many blocks as there are constraints. blocks is the partition of Asprion, Chinellato and Guzzella: a direct collocation transcription orders its variables by stage, so the Lagrangian Hessian is close to block diagonal in contiguous blocks and the block count is the stage count. Set partitioned_block_size to what one stage contributes (states × collocation points, plus controls) — too small and the block misses genuine intra-stage coupling, too large and each block carries more parameters than its one curvature pair per iteration can determine. The ordering assumption is reported rather than trusted: POUNCE_PARTITIONED_ORACLE prints the fraction of the exact Hessian’s Frobenius mass that falls inside the block pattern.

partitioned_update_type. SR1 is the default because an individual constraint is not convex. Damped BFGS would force every element model PSD, the solve would then scale it by a multiplier of either sign, and the indefiniteness would never reach the inertia correction. elements=blocks defaults to damped BFGS instead, since there the element is the Lagrangian restricted to a block, which is the object an interior-point method wants a positive-definite model of.

partitioned_max_element. An element with k nonzeros costs k(k+1)/2 stored reals, so one wide constraint row would dominate the memory. Elements wider than this degrade to a diagonal approximation satisfying the weak secant condition rather than being dropped, so a separable objective is still represented exactly and a coupled one approximately.

partitioned_curvature_cap — off, and every finite value measured was worse than off, non-monotonically so. On benchmarks/large_scale laptime at N=80 with max_iter=1200 (true optimum 65.462928): cap=1e1 exits ErrorInStepComputation at 1071 iterations and 65.518586; cap=1e2 hits the iteration limit at 67.202124; cap=1e6 hits it at 80.398129; off converges in 559 iterations at 65.462802. Rejecting an update is selective — it drops exactly the elements whose curvature is moving fastest and leaves those blocks stale while their neighbours update, and the resulting inconsistent Hessian costs more than a uniformly noisy but coherent one. Kept as a knob so the effect can be re-measured against a different element decomposition. Do not enable it without measuring.

Limited-memory Hessian (L-BFGS) initialization

Under hessian_approximation=limited-memory the Hessian model is B = σ I + V Vᵀ − U Uᵀ. The rank-2 corrections come from the curvature history; σ is the diagonal they are built on, and limited_memory_initialization chooses the formula for it.

The default matches Ipopt’s scalar1. Note that it changed in #677: every earlier release used scalar2 and ignored this option entirely — it was registered but never read, so setting it had no effect and no warning. The two differ by σ_scalar2/σ_scalar1 = (yᵀy·sᵀs)/(sᵀy)², which is ≥ 1 by Cauchy–Schwarz and grows without bound as the curvature pair becomes ill-conditioned, so on a badly scaled problem they are far apart. If you are reproducing results from an older POUNCE, set limited_memory_initialization scalar2.

history-max (#818)

Every upstream rule reads the newest curvature pair. history-max is POUNCE’s, and reads all of them: it applies the scalar1 formula to each pair in the history window and keeps the largest.

σ is the curvature the model assigns to every direction outside the span of the stored pairs — the rank-2 corrections say nothing there. sᵀy/sᵀs is a Rayleigh quotient of the true Hessian along one step, so when the problem’s curvature spans orders of magnitude it is an arbitrary sample of the spectrum. Land near the small end and B understates the curvature of every unexplored direction by up to cond(H); the step is longer than the truth by that factor and the line search has to claw it back. Over-stating σ costs an iteration; under-stating it costs a backtracking sweep and then feeds a tiny s back into the history. history-max takes the conservative reading. Nothing is lost on the directions the model does know: the last rank-2 update enforces B s_last = y_last whatever B0 was.

This is not the default, because it wins where the window cannot span the spectrum and loses where it can. On #818’s separable quadratic f(x) = Σ (sᵢxᵢ − 1)² with s = 10^linspace(0, 4, n) (iterations to converge, max_iter 2000, * = did not converge):

Reach for it when a limited-memory solve is stalling with a long backtracking sweep every iteration on a problem you know is badly conditioned, and the memory (limited_memory_max_history) is small relative to the number of variables. A running maximum over the whole solve rather than over the window was measured too and is worse than scalar1 at every size — it never comes back down, so a stiff early transient stays in B0 long after the iterate has left it.

recalc_y under L-BFGS

A quasi-Newton dual step is computed from an approximate Hessian, so an L-BFGS solve can settle a feasible primal and still fail to drive dual infeasibility to tolerance. recalc_y yes re-estimates the equality and inequality multipliers by least squares on every iteration whose constraint violation is below recalc_y_feas_tol (default 1e-6), side-stepping the approximation. Each firing costs one extra augmented-system solve.

Ipopt’s option text says this is used by default with a quasi-Newton Hessian. POUNCE does not enable it by default, because doing so regressed 7 of 57 fixtures on the L-BFGS leg — re-estimating y every iteration also overwrites Newton multipliers that were converging perfectly well. Reach for it when dual infeasibility oscillates without descending while the objective and the primal have already settled; that is the shape it fixes.

σ cannot be observed directly, but the symptom of a badly chosen one is recognisable: a search direction much larger than the problem’s scale, primal step sizes collapsing to 1e-3 or below, primal infeasibility that barely moves, and dual infeasibility climbing by orders of magnitude while the objective drifts.

Backtracking trial steps and alpha_red_factor_min

When the filter rejects a trial step the line search shortens it and tries again. Upstream reduces by a fixed factor, alpha *= alpha_red_factor (default 0.5), so reaching a step of size α takes log(1/α) trial points — each one a full objective evaluation.

That is cheap when the Hessian model is roughly right and ruinous when it is not. Under hessian_approximation=limited-memory the model’s scale can be wrong by orders of magnitude in any direction its curvature pairs do not span, so the acceptable step can be α ≈ 1e-6 and every iteration spends 19–20 trial points walking down to it (#818).

POUNCE therefore picks the next trial step by fitting the quadratic through the barrier objective’s value and slope at the current iterate and its value at the rejected trial, and jumping to that quadratic’s minimizer — the textbook safeguarded backtracking step. Two clamps bound it:

Acceptance is unchanged — this only decides which α is tried next, so no step it proposes can be accepted that the fixed sequence would have rejected. Set alpha_red_factor_min equal to alpha_red_factor to restore upstream’s fixed geometric sequence; an explicit value is honoured on both Hessian paths.

The default is off for the exact-Hessian path because a Newton step’s length is meaningful — the acceptable α is normally within a couple of halvings of 1, so interpolation buys nothing — and because enabling it there was measured to move 3 of the 156 fixture-legs in scripts/sweep-fixtures.sh — eigena2 27 → 31, issue_508_infeasible_gap_1em4 441 → 580 to the same certificate, and an objective digit on hs13_bigstart — with nothing on that leg improving in exchange.

It engages only once the fixed sequence has already spent six trial points. The interpolation treats a long line search, and it is only harmless where the line search is long: one that accepts in two or three trials never had the problem, and interpolating into it swaps a step length the filter was about to accept for a different one — a trajectory change bought for nothing. The threshold is not an option; it is the constant ALPHA_INTERP_MIN_TRIALS in crates/pounce-algorithm/src/line_search/backtracking.rs, which carries the sweep that chose it, every row against one baseline. Interpolating from the first trial instead moves 12 fixture-legs, three of them to a status the baseline did not have; gating at six moves four, one of them a gain (cresc4, Restoration_Failed/1323 → Solve_Succeeded/281) and none of them a loss.

Six rather than five because of two measurements. deb7 changes verdict at a gate of 5 (Error_In_Step_Computation → Restoration_Failed) and keeps it at 6; and one line search in the race_starts regression suite reaches exactly 5–6 trial points, so a gate of 5 interpolates into it and reroutes a whole multistart race, costing 32% more solver evaluations on one model there.

Where it loses, and what to do about it. Over a 32-cell sweep of the #818 model family (n ∈ {4, 8, 12, 20} × cond ∈ {1e2, 1e4, 1e8, 1e12} × limited_memory_max_history ∈ {6, 10}) the same 22 cells converge before and after — no cell gains or loses Solve_Succeeded — but it is not free: of those 22, 13 take fewer iterations, 5 are unchanged and 4 take more. The worst is n = 8 at cond 1e4 with memory 10, 188 → 580, a 3.1× regression to the same answer.

If a limited-memory solve stalls where it used to crawl, the first thing to try is limited_memory_max_history 10 — the cells this change does not fix are bounded by the quality of the quasi-Newton model, not by the trial sequence, and a wider window is what moves that bound. The 8-variable cond-1e8 case is the example: at memory 6 it exhausts max_iter — reaching f ≈ 9e-13 and x to 6e-7 relative, so it has found the answer but cannot certify it — and at memory 10 it converges in 61 iterations. The general escape hatch is alpha_red_factor_min equal to alpha_red_factor, which collapses the clamp and restores upstream’s sequence exactly.

One cell of that sweep, n = 8 at cond 1e12 with memory 6, read as a status regression during review — Diverging_Iterates at 352, at every alpha_red_factor_min and every gate measured. The cause was not the line search: the divergence guard was pronouncing unboundedness on a watchdog trial iterate, a point the line search had already rejected and was holding a snapshot to revert to. That is fixed separately. The cell still does not converge — it reports Error_In_Step_Computation at 521 iterations — but it no longer claims divergence, and it gets to a better objective than upstream’s sequence does in four times as many iterations (6.4e-11 against 2.8e-10 at max_iter).

When the line search fails anyway: limited_memory_ls_failure_restarts

When no trial step is acceptable, either the point is bad — infeasible, and the restoration phase is exactly the right tool — or the direction is, because W is a quasi-Newton model carrying curvature the iterate has left behind. Upstream has one answer for both, because restoration is the only fallback it has.

At an already-feasible point that answer is a no-op. The restoration NLP minimizes the constraint violation and there is none to minimize, so it wanders at θ ≈ 1e-13 and reports Restoration_Failed. On the deb7 fixture under limited-memory the solve stalls at inf_pr ≈ 1e-12 with inf_du ≈ 1e5, enters restoration at a point feasible to 8e-13, and spends 340 of its 1242 iterations there. On an unconstrained model θ is identically zero, so restoration cannot move at all.

The newest pair is kept rather than the history cleared, because σ is read off the history and an empty one falls back to limited_memory_init_val — a bare 1.0, which throws the model back to its first-iteration state on a problem whose curvature the solver has by now measured. This is L-BFGS-B’s col = 0 restart, adapted.

It is a rung and not a refusal: it fires only where restoration has nothing to reduce, it runs after the acceptable-point decline (so a point that already passes the acceptable tolerances is still reported rather than re-anchored and continued), and every path that reached restoration before still reaches it once the rung is spent. The bound is structural as well as counted — the re-anchor gives up once the history is down to one pair, so a second failure at the same iterate falls straight through. It has no effect under an exact Hessian, which has no curvature history to re-anchor.

It ships off, and it is not what fixes #818. The safeguarded interpolation above is; the rung was measured separately on top of it and does not pay for itself across the fixture corpus. Turning it on moves six lbfgs legs: deb7 715 → 610 iterations and issue_508_infeasible_gap_1em4 79 → 76 in its favour, against eigena2 91 → 98, pooling_rt2stp 295 → 307 and infeasible_square_scaled_1em4 24 → 26 — the last two being models the shipped configuration leaves exactly where main had them, so the rung introduces those two regressions rather than inheriting them. It stays in the tree, and stays documented, because the failure mode it treats is real and reproducible: on a model that stalls at inf_pr ≈ 1e-12 with inf_du large, limited_memory_ls_failure_restarts 1 is worth trying before concluding the solve is stuck.

Note that setting it — to any value, 0 included — opts the solve out of the automatic Solved_To_Acceptable_Level re-solve ladder, like every other option in TERMINATION_POLICY_OPTIONS. Leaving it unset is therefore not the same as passing limited_memory_ls_failure_restarts 0.

ℓ₁ penalty-barrier wrapper options

These tune the degenerate-NLP wrapper described in Running Solves. All are default-tuned and rarely need overriding:

The wrapper solves an augmented problem, c(x) − p + n = target with p, n ≥ 0, whose equality rows the slacks satisfy to machine precision by construction. So neither the residual the inner solve converged nor the slack sum Σ(p + n) is the violation of the constraints you declared — that is |pᵢ − nᵢ| per row, and at the barrier’s interior both slacks stay positive where their difference is zero.

Since gh#794 the wrapper therefore measures the original model’s rows and bounds at the returned point, and judges that violation by the tolerances you set — tol for a strict Solve_Succeeded, acceptable_tol for Solved_To_Acceptable_Level, scale-relative in both cases — before the ρ loop stops or the honest-infeasibility upgrade fires. l1_slack_tol survives only as the fallback for a model whose rows cannot be evaluated at the returned point, and Σ(p + n) keeps its other job as the Byrd-Nocedal-Waltz steering signal for ρ escalation.

That measurement is in your model’s own units, so it is reported in final_unscaled_constr_viol (and folded into final_unscaled_kkt_error), which is where the statistics contract puts original-unit residuals. The scaled family — final_constr_viol, final_kkt_error — carries the same number only when no row scaling is active, the case in which the two families are defined to agree anyway. Under an active nlp_scaling_method the scaled fields keep the augmented problem’s own residuals, because converting the measurement into that space needs row-scale factors the wrapper does not have once the inner solve has returned. The exit status is unaffected either way: it is decided from the measurement directly, not from these fields. If you are checking ℓ₁ feasibility programmatically, read final_unscaled_constr_viol.

NLP Presolve

POUNCE’s TNLP-wrapper presolve pipeline runs before the IPM starts. It tightens variable bounds, drops redundant rows, and (optionally) eliminates square auxiliary-equality sub-systems structurally. All are off by default — set the master switch first:

presolve=yes applies equally to CLI solves and to every IpoptApplication::optimize_tnlp library solve; callers no longer need to wrap a callback TNLP manually. The wrapper postsolves before finalize_solution, so callback payloads remain in the submitted TNLP’s original variable and constraint space. Bare callback TNLPs do not expose an expression provider. The pounce-rs builder can opt in through Problem::constraint_expression.

Linear-equality variable elimination (Phase 6)

presolve_linear_eq_reduction=yes is the only pass that removes columns. It reads the model’s linear equality rows and eliminates the variables they determine, iterating to a fixed point so chains propagate:

• a variable whose declared bounds are equal becomes a constant;
• a singleton row a·x = b pins its variable at b/a;
• a two-variable row a₁·x + a₂·y = b substitutes one variable for the other, x := α·y + β. There is no anchoring requirement: a row linking two otherwise-free interior variables — an arc equality, a Reference alias, a unit-conversion link — aggregates away, which is the case the auxiliary-equality pass cannot reach because it only solves determined square blocks.

Rows that collapse to 0 = 0 under the accumulated substitutions are dropped as structurally redundant.

A row written with a constant on the left — x0 − 2·x1 + 3 = 3 — is eligible on the same terms as x0 − 2·x1 = 0. The .nl reader folds a constant row body into the row’s bounds when the file is read, so the pass sees an ordinary linear equality.

Every eliminated variable’s bounds are transferred onto its survivor, so the reduced box is never looser than the original. finalize_solution lifts the primal back to the original variable order and recovers a multiplier for each consumed row, so .sol / JSON solution blocks keep the original model’s shape and can still be read positionally by AMPL or Pyomo.

Three things to know before turning it on:

• Dual attribution. A transferred bound’s multiplier comes back on the variable that declared the bound, not on the survivor that inherited it. The plan records where each reduced bound came from, and postsolve rescales the multiplier by the substitution’s coefficient α — and moves it to the other side of the box when α < 0, since a negative coefficient turns a lower bound into an upper one. On a model with a single active transferred bound the reported duals match a no-presolve solve exactly.

  Where the survivor’s own bound and a transferred bound are active at the same point, the split between the two multipliers is genuinely non-unique — the reduced problem has one where the full problem has two — and the pass leaves the whole multiplier on the survivor. That is a valid KKT point, but it is not the split a no-presolve solve happens to report. The same holds for a variable pinned by a singleton row a·x = b whose value lands on one of its own bounds: the row multiplier absorbs it.

  A practical consequence for .sol readers, unchanged by any of the above: the writer omits exact zeros from suffix blocks (it always has), so a variable whose bound multiplier is zero gets no ipopt_zL_out / ipopt_zU_out entry at all rather than an entry of zero. Code that indexes those suffixes must treat a missing index as zero — as it already must for any variable whose bound multiplier lands exactly on zero. Row multipliers are unaffected: the dual block is dense and comes back at the original row count.

• Bounds the reduced problem never saw. Re-attribution can only move a multiplier the solver reported, and sometimes there is none. The transfers can leave a survivor’s reduced box as a single point, and a variable with equal bounds is a fixed variable, which the solver drops — so it comes back with no bound multiplier at all even though the cluster it stands for is sitting on a bound that needs one. Postsolve fills that in: whatever stationarity residual the recovered row multipliers cannot close is a bound multiplier that was never reported, and it goes on the declared bound the point is actually resting on — the survivor’s own where that is the active one, otherwise the column the bound was borrowed from, through the same α rescale as above. A residual with no active declared bound to carry it is left alone rather than parked somewhere that would break complementarity.

  One column is deliberately outside this: one the model declares fixed (x_l == x_u). The solver drops those as parameters whether or not the reduction runs, and reports no multiplier for them either way, so nothing here changes what they report.

• Failing closed. If the equality system is contradictory, the pass stands down entirely and hands the model to the solver untouched, rather than being the first and only voice to call a model infeasible. The same goes for a model whose every column is determined: a zero-variable problem is not a shape worth handing the IPM.

It is off by default because it changes the variable count, which the sensitivity and reduced-Hessian paths index against the original .nl. (The CLI already disables presolve entirely when those are requested.)

LP and convex QP take a different route to the same reduction. Those models never reach Phase 6 — the CLI dispatches them to pounce-convex before any presolve wrapper is built — but they are not left unreduced. pounce-convex has its own presolve, on by default, and it now performs the two-variable aggregation as part of that catalog, sharing this planner rather than restating it. So the reduction is the same; only the switch differs (qp_presolve=no / presolve=no turns it off there, and presolve_linear_eq_reduction does not apply). See LP/QP routing.

The two agree on dual attribution as well: a transferred bound’s multiplier is reported on the column that declared the bound, not on the survivor that inherited it. They get there differently — Phase 6 records during planning where each reduced bound came from, while the convex path reads the leftover reduced cost at postsolve and hands it to whichever column is sitting on its own bound, because it also has inequality rows and its own bound-tightening layer to account for. Where the survivor’s own bound is active as well, both leave the multiplier on the survivor; the split is genuinely non-unique there and either answer is a valid KKT point.

Feasibility-based bound tightening (Phase 1b)

Interval-arithmetic propagation through nonlinear constraint expression DAGs (see FBBT). Available today for .nl-loaded problems via NlTnlp; other TNLP sources opt out silently.

Auxiliary-equality preprocessing (Phase 0)

A separate set of options controls the structural elimination pass documented in Auxiliary-Equality Preprocessing:

FERAL backend tuning

linear_solver=feral (the default — see Commonly used options) is configurable through thirteen feral_* options. Defaults are tuned for the IPM workload and rarely need changing; reach for these when profiling a specific problem. Each also falls back to a matching POUNCE_FERAL_* environment variable when left unset on the OptionsList (see Environment overrides).

feral_ordering variants

All six concrete and adaptive options live under the same string option. feral_ordering also falls back to the POUNCE_FERAL_ORDERING environment variable when not set on the OptionsList.

When in doubt: leave feral_ordering at the default. When a hard problem looks linear-solver-bound, try feral_ordering auto_race before per-variant manual sweeping — it’s the safe choice when the per-problem winner is uncertain.

Caller-supplied ordering (External)

Beyond the string variants above, a structure-aware caller can inject a precomputed permutation the generic AMD/METIS pass cannot see — a block-triangular / Schur ordering (Parker, Garcia & Bent, arXiv:2602.17968) or a tearing ordering from equation-oriented decomposition. Because a permutation is a vector it cannot travel through the string feral_ordering option; supply it programmatically instead:

• Python: Problem.set_ordering(perm) (and get_ordering() / clear_ordering()) — see the Python guide.
• Rust: IpoptApplication::set_external_ordering(perm).

perm is a 0-based, new-to-old permutation (perm[k] is the original index that becomes index k) whose length must equal the augmented KKT system dimension (variables + slacks + constraint duals), not the problem’s n. FERAL validates it as a bijection and fails the factorization with an error on a wrong length or duplicate — a valid but poor ordering only costs fill/time, never correctness. This maps to FERAL’s OrderingMethod::External (feral#107) and honors only the default FERAL backend.

Environment overrides (FERAL and debug gates)

A handful of knobs are reachable through environment variables. The feral_* numerics knobs read their POUNCE_FERAL_* variable only as a fallback when the matching option is left unset on the OptionsList — set the option (per solve, recordable, discoverable via the debugger’s opt command) in preference to the env var (process-wide, invisible to the solve report). The debug gates below have no option equivalent; they exist purely to switch on extra diagnostic output.

FERAL numerics fallbacks

Each maps one-to-one to a registered option in FERAL backend tuning. Prefer the option; the env var is the fallback for callers with no OptionsList (some tests, legacy embeddings).

These variables are parsed by feral::env (feral#176), which accepts the same spellings the option parser does — including scientific notation, so POUNCE_FERAL_MIN_PAR_FLOPS=1e8 sets the documented default rather than being silently discarded, as it was before feral 0.17.0 when pounce parsed these with a bare str::parse. A value that is out of range for the target type is clamped to the maximum; a value that cannot be parsed at all, or that fails the knob’s stated requirement (e.g. a negative pivot floor), is refused with a one-time warning on stderr and the default is used. A refused variable never silently changes numerics.

FERAL_PARALLEL (legacy, no POUNCE_ prefix) forces feral’s internal factor serial or parallel process-wide — 0/off/false/no to force serial, 1/on/true/yes to force parallel, and unset to leave feral’s own platform-derived default alone. The force-on direction is the only override available to CLI, Python and NL callers on a host where that autodetection is wrong (feral falls back to sequential when the rayon pool fails to build); the first-class per-backend lever, FeralConfig.parallel, is the Rust solver API, not an option.

Debug and diagnostic gates

These switch on extra diagnostic emission for a specific subsystem. Most emit at debug level under a pounce::* tracing target, so setting the gate alone is not enough — pair it with a matching RUST_LOG (e.g. RUST_LOG=pounce::mu=debug) or the output stays filtered. Presence-only unless a value is noted; they are diagnostic aids, not part of the stable interface, and may change between releases.

Two already-documented gates round out the set: POUNCE_DBG_LLM and POUNCE_DBG_VIEWER (see the debugger guide).

Logging and colored output

POUNCE emits structured logs and a colored iteration table through the tracing ecosystem. Behavior is governed by environment variables (not solver options), so they apply to the pounce CLI, the C/Python frontends, and anything embedding the library.

Filtering by subsystem. Solver internals log under namespaced targets — pounce::algorithm, pounce::linsol, pounce::mu, pounce::sqp, pounce::linesearch, pounce::restoration, pounce::presolve, pounce::py. For example, to trace only the restoration phase:

RUST_LOG=pounce::restoration=debug pounce problem.nl

Program output vs. logs. The iteration table, the final summary, and --dump diagnostics are program output on stdout; diagnostic and progress messages are logs on stderr. Redirecting one does not affect the other:

pounce problem.nl > result.txt 2> solve.log

Color. The iteration table is colored with a tiger/rust theme: restoration lines take a background that varies by restoration kind (soft-stay → tan, soft-exit → amber, hard → deep rust), and the row text shades from black toward red as the primal step length alpha shrinks (stalling). Color is emitted only when stdout is a terminal; redirected output and NO_COLOR get plain text with identical column alignment.

Machine-readable iterations. POUNCE_LOG_FORMAT=json turns the per-iteration records into JSON on stderr:

POUNCE_LOG_FORMAT=json pounce problem.nl 2> iters.jsonl

LP / QP Solver Routing

POUNCE can route linear programs (LP), convex quadratic programs (QP), and convex quadratically-constrained QPs (QCQP) to a specialized interior-point solver (pounce-convex) instead of the general nonlinear (NLP) filter-IPM. The specialized path uses Mehrotra predictor-corrector and reaches the solution in materially fewer iterations on these problem classes — typically 30–50% fewer than the general NLP path on bound- or inequality-constrained convex QPs.

Routing is automatic and transparent: you do not change how you call POUNCE. The same pounce problem.nl, the same SolverFactory('pounce') in Pyomo, and the same AMPL solve all work unchanged — POUNCE inspects the problem and picks the solver.

How routing works

When POUNCE loads a problem it classifies it into one of:

The classifier is conservative: a problem is sent to the convex solver only when POUNCE can prove it is convex — an LP or convex QP (degree-≤2 objective with a positive-semidefinite Hessian, linear constraints), or a convex QCQP (additionally allowing convex-quadratic inequality constraints, each with a positive-semidefinite Hessian and a one-sided ≤ bound, which are reformulated to second-order cones). Anything it cannot prove convex — transcendental terms, an indefinite objective Hessian, a quadratic equality, or a quadratic inequality whose feasible set is nonconvex — falls back to the general NLP solver, which always produces a correct (locally optimal) answer. You never get a wrong “optimum” from a misclassification.

The nonconvex QP class is narrower than “any quadratic that is not convex”: it means an indefinite objective Hessian over linear rows. A model that is also curved in its constraints is a nonconvex QCQP and classifies NLP, because the QP extractor behind the convex path keeps only the linear part of each row — calling such a model a QP would hand an engine a problem with its curved constraints deleted. POUNCE_DBG_CLASSIFY=1 names which of the two a model landed in.

Note on QP detection. The AMPL .nl format has no dedicated quadratic section: a QP’s quadratic terms are written into the nonlinear expression tree. POUNCE walks that tree to recover the Hessian and test convexity, the same way QP-capable AMPL solvers do.

Note on row constants. A .nl writer may leave a constant on the left of a constraint — x0 + x1 + 3 <= 6 rather than x0 + x1 <= 3 — and it too lands in the nonlinear expression tree. The reader folds such a constant into the row’s bounds when the file is read, so a model that is otherwise an LP still classifies as one. The shift is exact: body and bound move together, so the solution and every multiplier are the same as for the hand-folded model.

Choosing the solver explicitly

The solver_selection option overrides the automatic choice. It is a normal POUNCE option, so it works on the command line, in an options file, or through Pyomo’s solver.options.

# Let POUNCE decide (default):
pounce model.nl

# Force the NLP path even on a convex QP (e.g. to compare):
pounce model.nl solver_selection=nlp

# Insist the problem is a convex QP — fail loudly if it is not:
pounce model.nl solver_selection=qp-ipm

# Solve that same QP with the active-set engine instead of the IPM:
pounce model.nl solver_selection=qp-active-set

A forced value that does not match the detected class is rejected with a clear message rather than silently ignored:

pounce: problem class NLP does not match forced solver qp-ipm
        (expected an LP or convex QP)

qp-active-set hands the QP directly to pounce-qp’s ParametricActiveSetSolver, through the same convex driver the IPM uses — so it inherits presolve, postsolve, dual recovery, .sol writing, timing and the convex status vocabulary. It is not the same route as algorithm=active-set-sqp, which wraps the QP in the full SQP outer loop; that option still exists and is the right one for a genuine NLP.

It is the one forced selection that takes a nonconvex problem. pounce-qp handles an indefinite Hessian by construction — inertia control shifts the H block until the reduced KKT factor has the right inertia — so solver_selection=qp-active-set on a nonconvex QP solves it rather than refusing it. What comes back is then a local optimum, exactly as Optimal Solution Found means locally optimal on the NLP path, and the status line names the class so the reader can tell:

Problem class: nonconvex QP. Selected solver: active-set QP (pounce-qp) [solver_selection=qp-active-set].
POUNCE (nonconvex QP active-set, pounce-qp): Optimal Solution Found.  obj=...

auto still sends that class to the NLP filter-IPM. The class is POUNCE’s inference, and for a nonconvex model the general path is the safer default — so the active-set route is reachable only by naming the engine. A nonconvex QCQP is refused on this route for the reason given under How routing works above: it does not classify as a QP at all.

Choose it deliberately. For a cold, one-shot convex QP the interior-point path (qp-ipm, and what auto selects) is materially more robust: on the 138-problem Maros-Mészáros set the IPM solves 137 while the active-set engine solves substantially fewer, mostly by exhausting its iteration budget on large degenerate instances. That is the expected character of a cold active-set method rather than a defect — its iteration count is combinatorial in the size of the active set, where an interior-point count is nearly independent of problem size. The active-set engine earns its keep on warm-started sequences — MPC steps, branch-and-bound nodes, continuation — where consecutive QPs differ little and the working set carries over; see solve_parametric.

What it will not do is lie: it reports Maximum_Iterations_Exceeded rather than a wrong answer, and every claimed optimum is re-verified against the original problem’s KKT conditions before being reported.

From Pyomo

solver = SolverFactory('pounce')
solver.options['solver_selection'] = 'qp-ipm'   # or 'auto', 'nlp', ...
solver.solve(model)

What you get back

Before solving, POUNCE prints a one-line routing banner naming the detected class, the solver it selected, and the effective solver_selection — so it is always clear which of POUNCE’s solvers ran and why:

Problem class: LP. Selected solver: convex QP interior-point (pounce-convex) [solver_selection=auto].

(The banner is suppressed alongside the startup banner — sb yes or JSON-debug protocol mode — to keep stdout clean for machine consumers.)

The convex IPM then reports the same way as the NLP path: an optimal-status line, the objective value (in your original sense — a maximize objective and any constant term are reported correctly), and a .sol file with the primal solution when one is requested.

POUNCE (LP IPM, pounce-convex): Optimal Solution Found.
        obj=2.00000000  iters=2

Driver. The convex path uses the homogeneous self-dual embedding (HSDE) interior-point driver — the same self-dual formulation Clarabel/ECOS use. It is self-starting, returns verified infeasibility/unboundedness certificates, and conditions the KKT system internally through its per-cone scaling, so it solves even badly-scaled LPs (e.g. NETLIB nl, ‖c‖ ~ 1e6) without external pre-scaling.

Presolve

Before the convex interior-point solve, POUNCE runs a presolve pass that shrinks the problem and can detect trivial infeasibility or unboundedness without solving. It removes empty, duplicate, and activity-redundant rows; fixes and substitutes structural columns (singleton-row fixings, free columns, free column singletons); folds away two-variable equality rows (below); and recovers both the primal and dual of the eliminated pieces so the reported solution is for your original problem. When it reduces the model, it logs a one-line summary:

Presolve: 40 → 24 vars, 12 → 4 rows (fixed 3, free-fixed 2, substituted 3, aggregated 8, ...)

Two-variable equality rows (aggregation)

A row a₁·x + a₂·y = b linking two variables says one of them is the other, up to a scale and a shift — an arc equality between two units, a Reference alias, a unit conversion. Neither variable is determined by it, so nothing in the older catalog could act on it, and on a flowsheet these rows are most of the model. POUNCE now substitutes one variable for the other and drops the row, iterating to a fixed point so chains of aliases collapse to a single column. Any bound on the eliminated variable is carried across onto the one that survives, so the reduced problem describes exactly the same feasible set.

Two things this deliberately does not do:

• It never calls your model infeasible. A contradictory alias system — x = y and x = y + 1 — makes the pass stand down and hand the model over untouched, for the rest of presolve or the solver itself to judge.
• It does not run on the conic path (SOCP, exponential/power cones, SDP, SOS). Those rows are structurally coupled in fixed-size blocks that a substitution would rewrite.

The aggregation shares its planner with the NLP path’s Phase 6, so the two agree on what can be eliminated (see NLP Presolve).

Infeasibility verdicts are re-derived before they are reported

A presolve infeasibility comes back in milliseconds with no iteration behind it, so when it is wrong it is the most expensive answer the solver can give. Two reductions — forcing constraints and dominated columns — fix a variable at a value they choose from a tolerance judgment, and a fixing that is wrong is substituted into every row that variable appears in until some row reads as contradictory: a false infeasibility, reported against a row nowhere near the reduction that caused it.

So presolve does not report an infeasibility on the strength of the pass that found it. It re-derives the verdict from your original model with those two reductions switched off, and reports Infeasible_Problem_Detected only if that pass reaches the same conclusion on its own. If it does not, the model is solved normally and presolve says so:

Presolve: discarded an unconfirmed infeasibility claim — <screen> (<detail>); solving normally

A confirmed verdict now names the screen that proved it and the row, column, or bound it tripped on, rather than exiting silently:

Presolve: proved primal infeasible — empty equality row (equality row 7 is `0 = 3e0`)

Nothing that only reports is withheld from the re-derivation — empty rows, activity ranges, parallel rows, and emptied-row residuals all still apply — so no infeasibility presolve could detect before goes undetected now. What the guard costs, in the rare case it fires, is a handful of eliminations.

When the reduction is truncated

The reductions are iterated to a fixpoint — each one can expose work for the next, so presolve keeps going until nothing fires. It also carries a cap on how many layers that may take, and on a model with a long bound-propagation chain the cap is what stops it. When that happens the summary line says so:

Presolve: 315 → 128 vars, 233 → 77 rows (fixed 61, ..., tightened 158, cap-truncated after 32 layers)

This is common and it is not a problem. Measured across the LP and QP suites, the cap binds on 46% of LP models and 25% of QP models — and on every one of the 394 models that presolve at all, it changed only how tightly variable boxes were narrowed, never the structural reduction: same variables, same rows, same fixings, aggregations, forcing rows and dominated columns as running the iteration to convergence. Bound propagation is the one reduction that can keep going indefinitely, so it is what the cap ends up trimming.

What you get is still a correct problem — every reduction applied is a sound transform with its own dual recovery, and your solution is postsolved back to the original either way. The suffix is there so a reduction that came out of a truncated loop is distinguishable from one that converged, which matters when you are comparing two runs or reporting a bug against presolve. There is no option to turn it up.

Presolve is on by default. Turn it off with qp_presolve=no (e.g. to compare timings or isolate a solver issue):

pounce model.nl qp_presolve=no

Presolve on a convex QCQP

The switch applies to the conic driver too — the one that solves convex QCQPs. Two things about it differ from the LP/QP path, and both follow from the same fact: a quadratic constraint is reformulated into a second-order cone block, and a cone block’s rows are coupled to one another.

Only the ordinary linear rows are reduced. Every row of a cone block is protected: it is never dropped, never merged with another row, never used to tighten a variable bound, and the variables it couples are excluded from the dominated-column reduction. Dropping any single row of a block would change which constraint the block encodes, with nothing to signal it — the answer would simply come back wrong. So on a model that is a variable box plus quadratic constraints and nothing else, presolve has nothing to act on and prints no summary line. That is the expected result, not a failure.

The loop runs once, not to a fixpoint. The reduced cone partition has to be readable off the surviving rows, which holds for a single pass. There is no cap-truncated suffix on this path for the same reason.

Where a QCQP does carry ordinary linear inequalities — which is the common shape — those are reduced exactly as on the LP/QP path, and the summary line looks the same:

Presolve: 4 → 4 vars, 9 → 7 rows (fixed 0, free-fixed 0, substituted 0, forcing 0, dominated 0, tightened 0)

Tuning the convex IPM

Beyond the shared tol and max_iter, the convex engine takes these:

qp_gondzio_corr is worth a sentence on where it does and does not apply. The correctors box-project the complementarity products sᵢzᵢ back into [0.1·μ, 10·μ], which needs the product to be elementwise — so the loop is gated on the cone being a pure nonnegative orthant and a solve carrying a single second-order or PSD block never enters it. That includes convex QCQP on the conic route, whose whole point is the SOC reformulation. POUNCE_DBG_GONDZIO=1 prints one line per convex solve — iterations, correctors attempted, correctors accepted and the mean step gain — which is the direct way to check whether the scheme is doing anything on your model before tuning the number.

These reach the engine through the pounce CLI, which is the one entry point that classifies a .nl model and routes it. A library solve refuses a non-default value rather than accepting one it would drop — IpoptApplication has no structure extraction, so it cannot route to the convex engines at all (the same reason solver_selection=lp-ipm errors there). From Python, pounce.solve_qp / pounce.solve_cone drive the engine directly and take these knobs as typed arguments.

Scope and limitations

• Convex problems only. Nonconvex (indefinite-Hessian) QPs, quadratic equalities, and quadratic inequalities whose feasible set is nonconvex are solved by the NLP path to a local minimum; POUNCE does not do global optimization.
• Convex QCQP (convex-quadratic constraints) routes to the conic IPM: each convex-quadratic inequality ½xᵀQx + aᵀx + b ≤ 0 (with Q ⪰ 0) is reformulated to one second-order cone (Q = FᵀF, so ‖Fx‖² = xᵀQx) and solved alongside the QP objective and linear constraints.

Both the primal solution and the constraint duals are written to the .sol file, in the same sign convention as POUNCE’s NLP path (so Pyomo and AMPL read them identically regardless of which solver ran).

Requests the convex path does not implement

The convex solvers are a specialized fast path, not a drop-in for every option the NLP path honors. Where a request would be dropped rather than merely unused, routing gives way rather than answering a different question:

A positive obj_scaling_factor is not in this table: it only rescales conditioning, and the convex path reports natural units either way, so both paths give the same answer.

Sensitivity is served here, not routed away

A plain parametric sensitivity request — the sIPOPT sens_* suffixes, without a reduced-Hessian request — is no longer in the table above. On an LP or convex QP the convex path computes it directly, through pounce_convex::QpSensitivity, and writes the same sens_sol_state_1 block the NLP path writes; a .sol consumer cannot tell which engine answered, and the banner names the one that did.

This was a reroute until the convex arm grew a parametric step of its own. It is worth knowing which way it goes, because the two engines have very different costs on a large LP and the answer is the same either way.

Three things keep it in the table:

• A reduced-Hessian request. QpSensitivity has one, but it is a different computation behind the same word — a null-space projection where the CLI’s sIPOPT path takes the Schur route. Serving it here would silently change which number --compute-red-hessian returns.
• A conic model. Every cone family has a face decomposition now (see The convex/conic solver), but the CLI’s conic dispatch extracts through its own provenance map and mapping pins through that is unwritten.
• A pin that is not x_p = p₀ with a unit coefficient. The convex step perturbs the equality right-hand side b; an inequality pin lives in h, which is a different perturbation. Rather than answer a different question, the model goes to the path that has always handled it.

Presolve is switched off for a run that serves a sensitivity request, on this path exactly as on the NLP one — but not for the reason it looks like. The convex driver postsolves back to the extracted-QP space before anything downstream runs, so the pins stay valid with presolve on and the step is still within 1e-6 of the NLP path’s. What presolve costs is accuracy: it can fix the very parameter the pin parametrizes and drop its row, leaving the sensitivity to read a postsolve reconstruction rather than the KKT the solve converged — four orders on the fixture that exercises it. Whether that can also move the active set the sensitivity infers, which would be a wrong derivative rather than a less accurate one, is not yet measured.

When the convex path cannot certify an LP

Routing gives way one more time, and this one is decided after the solve rather than before it. Under auto, an LP whose convex solve finishes without a KKT certificate — Solved to acceptable level (reduced accuracy) or Maximum iterations exceeded — is re-solved on the general NLP interior-point path, which owns the whole verdict. Nothing from the declined convex solve is printed or written, so a rerouted run still reports exactly one status.

The case this exists for is the NETLIB gen / gen1 family. They are highly degenerate and rank-deficient, strict complementarity fails, and a pure interior-point method cannot certify the optimal vertex: the convex IPM spends its whole 200-iteration budget (190.8 s) and stops at a primal residual of 1.4e-7 against tol = 1e-8. The NLP filter-IPM — the same binary, the default for every other class — solves the same model in 19 iterations and 0.98 s to a strict certificate, matching Ipopt-3.14.20/MA57 to four figures. Rerouting is also the faster answer here: a second solve of one second is nothing against the three minutes the first one costs.

The fallback is narrow by construction, and does not fire when:

A tightened tol is deliberately not in that list: that is an accuracy request, so trying the engine that can meet it is the right response.

An explicitly set max_wall_time is forwarded to every automatic convex LP, QP, active-set QP, and SOCP route. Time spent extracting and presolving the convex model is charged to the same budget as all engine retries. Expiration is reported as MaximumWallTimeExceeded, AMPL solve_result_num = 400, with the message “Maximum wallclock time exceeded.” A timed-out convex solve is final: automatic routing never starts a fresh convex attempt or falls back to the NLP engine. max_cpu_time remains an NLP-side option and is not forwarded.

Infeasible and unbounded problems

The convex solver detects infeasibility and unboundedness directly, reporting a clean status instead of exhausting the iteration budget:

• Primal infeasible — no point satisfies the constraints. Reported with AMPL solve_result_num 200.
• Unbounded (dual infeasible) — the objective decreases without bound along a feasible direction. Reported with solve_result_num 300.

Each verdict is backed by a verified certificate (a Farkas infeasibility proof or an unbounded recession direction that is checked, not merely inferred), so these statuses are never reported in error; a problem the solver cannot certify simply runs to the iteration limit — and, if it is an LP under auto, is then handed to the NLP path (above).

solver_selection=qp-active-set follows the same contract. Its inner QP certifies the recession ray of the linearization, which on a nonlinear model is not yet a statement about the problem, so the ray is re-tested against the true objective and constraints before the 300 is reported; a ray that does not survive yields Search_Direction_Becomes_Too_Small, never an unboundedness claim.

The design and roadmap live in dev-notes/lp-qp-routing.md.

Convex Solver: LP, QP, and SOCP

POUNCE ships a specialized convex conic interior-point solver (pounce-convex) alongside the general NLP filter-IPM. It solves the standard-form convex program

minimize    ½ xᵀP x + cᵀx
subject to  A x = b
            G x ⪯_K h
            lb ≤ x ≤ ub

where P ⪰ 0 and the inequality block lies in a product cone K of nonnegative orthants and second-order cones. P = 0 is an LP; an all-orthant K is an LP/QP; second-order blocks make it an SOCP.

The method is a Mehrotra predictor–corrector primal–dual interior-point algorithm with Nesterov–Todd scaling for the cones, sharing the pure-Rust feral sparse LDLᵀ backend with the NLP path. It reaches optimality in materially fewer iterations than routing the same problem through the general NLP solver (≈30–50% fewer on bound/inequality QPs).

Inspiration. The conic interior-point design follows Clarabel (Goulart & Chen) — handling a quadratic objective directly and a product of symmetric cones — and the presolve follows PaPILO (the presolving library of SCIP). POUNCE does not wrap either (the pure-Rust guarantee) but ports their ideas; see Acknowledgments.

This chapter covers the Python API (pounce.qp and the differentiable pounce.jax layers). For automatic CLI/Pyomo routing of .nl LPs/QPs, see LP / QP Solver Routing. Runnable, progressive notebooks live in python/notebooks/: 15_convex_qp.ipynb, 16_socp.ipynb, 17_differentiable_convex.ipynb.

Quadratic programs

import numpy as np
from pounce.qp import solve_qp

# min ½·2‖x‖² − 3x₀ − 4x₁  s.t.  x₀ + x₁ ≤ 1,  0 ≤ x ≤ 1
r = solve_qp(
    P=np.diag([2.0, 2.0]),
    c=[-3.0, -4.0],
    G=[[1.0, 1.0]], h=[1.0],
    lb=[0, 0], ub=[1, 1],
)
r.status   # 'optimal'
r.x        # primal solution
r.y, r.z   # equality / inequality multipliers
r.z_lb, r.z_ub  # bound multipliers (≥ 0)
r.obj, r.iters

P (lower triangle used, assumed symmetric), A, and G accept dense arrays or scipy-sparse matrices; any of them may be omitted. The result is a QpResult dataclass with a .success property. The solver reports verified infeasibility / unboundedness ('primal_infeasible' / 'dual_infeasible') backed by a Farkas / recession certificate rather than an iteration-limit guess.

Second-order cone programs

A second-order (Lorentz) cone is { (t, x) : t ≥ ‖x‖₂ }. Partition the inequality rows of Gx ⪯_K h with cones — a list of (kind, dim) specs ("nonneg" or "soc"; a bare int means a second-order cone). Each slack block s = h − Gx must lie in its cone.

from pounce.qp import solve_socp

# minimize ‖x − x*‖  ⇔  min t s.t. (t, x − x*) ∈ SOC
r = solve_socp(
    c=[1.0, 0.0, 0.0],                 # minimize t
    G=-np.eye(3), h=[0.0, -2.0, 1.0],  # s = (t, x₀−2, x₁+1) ∈ SOC(3)
    cones=[("soc", 3)],
)
r.x   # ≈ [0, 2, -1]:  t* = 0, x = x*

Mixed cones compose — e.g. cones=[("nonneg", 1), ("soc", 2)] puts the first slack in ℝ₊ and the next two in a 2-D second-order cone. Large cones use a sparse diagonal-plus-rank-1 KKT representation (one auxiliary variable per cone, the ECOS/Clarabel “sparse SOC” trick) so the factorization stays sparse.

Warm starting

Feed a previous (or nearby) solution back to seed the interior-point iteration — useful for parametric sweeps, receding-horizon MPC, and branch-and-bound subproblems:

base = solve_qp(P=P, c=c, G=G, h=h, lb=lb, ub=ub)
nxt  = solve_qp(P=P, c=c2, G=G, h=h, lb=lb, ub=ub, warm_start=base)

The warm start only affects the iteration count, never the solution (a mismatch is ignored). The recentering is adaptive for the orthant (sized to the warm point’s KKT residual, so it exploits a nearby problem’s duals yet self-corrects when the active set moves) and re-centers the cone duals for second-order blocks (a converged conic point sits on the cone boundary, where the scaling is singular).

The step length is what makes it pay off

A warm start lowers the starting duality measure μ₀; whether that turns into fewer iterations depends on how much of each Newton step the solver is allowed to take. With a static fraction-to-boundary parameter τ, every step covers at most a τ fraction of the distance to the cone boundary, so μ falls by a fixed factor per iteration and the count is log₁/₍₁₋τ₎(μ₀/tol) however good the start was — a logarithm of the perturbation, not the one or two Newton steps a nearby problem deserves.

So on orthant blocks the step follows the Mehrotra tail τ = clamp(1 − μ, tau, tau_max): as the solve converges τ approaches 1 and a near-optimal iterate takes a near-full Newton step. On the QP families in the warm-start benchmark this is worth 35–60% of the warm iterations. Both ends are tunable, and both are method="ipm" only:

r = solve_qp(P=P, c=c2, G=G, h=h, warm_start=base,
             tau=0.95,        # floor: the flat τ far from the solution
             tau_max=0.999)   # ceiling on the tail (default: just under 1)

Passing tau_max=tau pins τ flat — the most conservative setting, and the one to reach for if a badly-conditioned sequence starts producing numerical_failure. Two scopes are deliberate and not tunable: second-order and PSD blocks always keep the static tau (their boundary is curved, and an iterate that close to it breaks the Nesterov–Todd scaling), and cold solves are unaffected because they run the homogeneous self-dual embedding, a different loop.

Wall-clock budgets

solve_qp, solve_socp, solve_qp_batch, and solve_qp_multi_rhs take a time_limit in seconds (None, the default, means unbounded). Reach for it when an answer is needed on a schedule — a receding-horizon controller with a fixed control period, a sweep where one pathological instance must not stall the rest, or any solve sitting behind a request:

r = solve_qp(P=P, c=c, G=G, h=h, warm_start=previous, time_limit=0.005)
if r.status == "time_limit":
    ...   # `r.x` is the best iterate reached, not a KKT point

max_iter cannot express this. One interior-point iteration may be a single KKT solve, or a factorization plus several inertia-controlled refactorizations with escalating shifts, and the LP route can add a simplex crossover phase — so per-iteration cost varies by more than an order of magnitude within one solve, before problem size enters into it. No iteration count means “5 ms” across two problems.

Three properties are worth knowing:

• A verdict outranks the clock. optimal, optimal_inaccurate, primal_infeasible, and dual_infeasible survive a deadline that passed while the solve was finishing; only a give-up result is relabelled time_limit. So the status is always truthful about what was proved, and a budget can never turn into a wrong optimal.
• The budget is per solve, not per call. On the batched entry points each instance opens its own deadline scope, so time_limit=10 over 100 problems permits 1000 s of wall clock. A shared clock would make which instances get cancelled depend on rayon’s scheduling, and so on the machine. Bound the whole call around the call.
• Results become machine- and load-dependent, inherently — which is why this is opt-in and absent from the default path. An in-flight factorization is not interrupted, so expiry can overshoot by one such operation.

The differentiable layers (pounce.jax, pounce.torch) deliberately do not take one: they raise on time_limit because a non-KKT iterate makes the implicit-function gradient meaningless, and silently wrong gradients under load are worse than a slow layer. On the CLI the same mechanism is spelled max_wall_time.

Batching and factorization reuse

from pounce.qp import solve_qp_batch, QpFactorization

# Solve many independent QPs in parallel (rayon, across instances).
results = solve_qp_batch([dict(P=P, c=c_k, G=G, h=h) for c_k in cs])

# Build the KKT symbolic factor once, solve many same-structure problems.
fac = QpFactorization(P=P, c=c0, G=G, h=h, lb=lb, ub=ub)
for c_k in cs:
    rk = fac.solve(P=P, c=c_k, G=G, h=h, lb=lb, ub=ub)  # reuses the factor

solve_qp_batch parallelizes across instances (outer-parallel / inner-serial) and QpFactorization reuses the AMD ordering and symbolic factorization across solves that share a structure — the two compose with warm starting.

Post-optimal sensitivity (QpSensitivity)

QpSensitivity is the convex arm’s sIPOPT analog: it holds the factored active-set KKT system at the optimum, so each parametric_step is a single back-substitution.

from pounce.qp import QpSensitivity

# min ½‖x‖²  s.t.  x₀ + x₁ = 2   →   x* = (1, 1),  dx/db = (½, ½)
s = QpSensitivity(P=np.eye(2), c=[0.0, 0.0], A=[[1.0, 1.0]], b=[2.0])
dx = s.parametric_step([0], [1.0])       # perturb b₀ by +1

It perturbs the equality right-hand side b, and reports the active set, the weakly-active set, a reduced Hessian, and two conditioning diagnostics (ill_conditioned, last_step_residual) that let a caller detect a step it should not trust.

Holding the step inside the bounds

The plain step is a linear predictor, so a large enough perturbation can point outside the variable box. parametric_step_bounded repairs that the way the NLP arm does — by pinning the crossing coordinate at its bound and re-solving, so the other coordinates move to suit and the constraints still hold. Clipping instead would satisfy the bounds and quietly break the equalities.

#![allow(unused)]
fn main() {
let (dx, pinned, stop) =
    sens.parametric_step_bounded(&[0], &[-6.0], /* bound_eps */ 1e-3, /* max_iter */ 16)?;
}

This is not a second implementation: it runs pounce_sens_core::boundcheck::refine_step_onto_bounds, the same code the NLP arm runs, reached through QpSensitivity::backsolver(). That machinery is generic over the SensBacksolver trait, whose whole required surface is dim() and solve(rhs, lhs), so an engine that can back-solve against its converged factor gets fix-relax, path following and the directional derivative without porting any of them.

Both halves of fix-relax are available: a coordinate the step carries past a bound is pinned there, and a bound whose multiplier the step drives negative is released so the variable can leave. Releasing is exact on this arm — the convex active-set KKT has no barrier term to destroy, so it costs one numeric refactorization against an unchanged sparsity pattern.

Following the path

parametric_step_path applies the perturbation a little at a time, stopping wherever the active set changes:

#![allow(unused)]
fn main() {
let (dx, segments) = sens.parametric_step_path(&[0], &[3.0], /* max_iter */ 32)?;
for s in &segments {
    println!("at {:.3}: x{} {} its {} bound",
             s.at, s.var_row, if s.pinned { "reached" } else { "left" },
             if s.lower { "lower" } else { "upper" });
}
}

A QP’s solution path is piecewise affine, so within a segment the walk is exact and the reported breakpoints are the real ones. Use this when a perturbation is large enough to change the active set more than once, or when you want the events rather than only the endpoint.

What each bound is doing

activity() classifies every bounded variable and every inequality row:

#![allow(unused)]
fn main() {
let rep = sens.activity();
rep.var_status[j]   // INACTIVE / WEAKLY_ACTIVE / STRONGLY_ACTIVE / AMBIGUOUS / …
rep.var_ratio[j]    // the ratio Σ/q the verdict came from
rep.mu              // the achieved complementarity it banded against
}

The rule is the same one the NLP arm applies — pounce_sens_core’s activity kernel — so the two arms agree on what a kink is. A kink is a bound whose slack and multiplier vanish together: there the derivative is two-valued and a parametric step needs a side, which is what the fix-relax and path modes above are for.

AMBIGUOUS is not “probably not a kink”. A genuine kink lands there whenever its coordinate is coupled to another through the Hessian, because the curvature the classifier can afford is a diagonal (for a variable) or the curvature along the row’s own gradient (for a row), while the multiplier is generated by the curvature reduced along that coordinate. The ratio is reduced/diagonal, which equals one only when the coordinate is decoupled — and it does not depend on μ, so solving more tightly will not resolve it. Never read the activity class as a proxy for kink-ness.

Degenerate LPs need crossover

lp_without_crossover() is true when the problem is a pure LP (P = 0) whose solve did not run crossover. At a degenerate optimal vertex more constraints are active than there are variables, the active-set KKT is rank-deficient, and dx/db is not single-valued — on a two-variable example the step comes back summing to half the perturbation it should. ill_conditioned() already catches that; this flag names the cause. The fix is to solve with qp_crossover=yes, so the interior point is pivoted to an exact vertex basis first.

Because the flag reads opts.crossover, the options you hand to build must be the options the solve actually ran with.

Orthant rows, and second-order cones

QpSensitivity::build covers LP and convex QP — problems whose inequality block is a nonnegative orthant. Cones go through build_conic instead, which handles SecondOrder blocks and refuses the others.

That distinction matters more than it looks, because solve_socp_ipm and solve_qp_ipm return the same QpSolution type and the cone partition travels beside it as a separate cones argument. So on the Rust API, handing a solved conic program to QpSensitivity::build used to be accepted and answered — every cone row read as an orthant row, producing a number that was not a derivative, with no warning. It is now refused with SensError::NotOrthantComplementary: an orthant row complements row by row (sᵢ ≥ 0, zᵢ ≥ 0, sᵢzᵢ ≈ μ), while a cone satisfies only the block inner product ⟨s, z⟩ = 0.

Python callers were never exposed to this: pounce.qp.QpSensitivity solves internally with the QP interior-point solver and accepts no cones=.

Cones: the face, not the rows

Use QpSensitivity::build_conic(prob, cones, sol, opts, active_tol, backend) for a problem that carries cones. An all-Nonneg partition is the orthant problem and delegates to build, so the two entry points cannot answer differently on the same input.

A cone’s active object is not a set of rows. Its slack sits on a face, and every family splits the same three ways — reported by cone_block_kinds() as ConeBlockKind:

What the boundary face is, and how many rows it contributes, is per family:

The PSD case is the one that is not just another smooth facet: its face has codimension q(q+1)/2, so a Psd(3) block at rank 1 contributes three rows. Its tangent is Vᵀ dX V = 0 for V a basis of ker S, which is the first-order form of the Schur complement C − Bᵀ A⁻¹ B vanishing.

There is no “unsupported cone” error. The match that dispatches the face decomposition is exhaustive over ConeSpec, so a family added later is a compile error rather than a runtime refusal — a stronger promise than a message, and it keeps an empty error category from sitting in the public API looking like a live one. What gets refused is a point, not a family.

The boundary curvature is part of the answer

Every orthant row and every variable bound is a hyperplane, so the sensitivity KKT’s (x,x) block is the objective’s Hessian P and nothing else. Every conic boundary face is curved, and its curvature enters the same block:

  second-order:   H = P + (ν/s₀) · ( Σ_{r≥1} gᵣgᵣᵀ − u uᵀ ),  u = Σ_{r≥1} (sᵣ/s₀) gᵣ
  exp / power:    H = P − ν · Gᵀ ∇²φ G                        (rank one, both)
  PSD at rank r:  H = P + 2 · Σ_{l ≤ r} Σ_{k ≤ q} (λ_k / a_l) · c_lk c_lkᵀ,
                  c_lk = Gᵀ svec(sym(ũ_l w̃_kᵀ))

with ν the multiplier on the facet’s φ, and for the PSD case a_l, ũ_l the slack’s positive eigenpairs and λ_k, w̃_k the dual’s. Every one of these is positive semidefinite, as a concave constraint’s contribution must be. This is not a refinement. Omit it and the step converges to the wrong derivative: on the worked fixture in crates/pounce-convex/tests/convex_soc_sensitivity.rs, dx/db reads (0.348, 0.652) where the closed-form answer is (0.5, 0.5), at every perturbation size, while every internal residual stays happy — the step solves exactly the KKT it was given, and that KKT is not the problem’s. The guard that catches it is the re-solve oracle in that file, the one test in the crate that compares against a number the sensitivity layer did not produce.

Where it refuses

Two errors mark the refusals, and the split between them is load-bearing. SensError::NonsmoothConePoint { block, what } means no single dx/db exists here — a kink, a collapsed normal, a two-valued derivative. SensError::ActiveSetOverdetermined { block, what } means the derivative exists and this active set cannot express it; a caller matching the first to decide “genuinely nondifferentiable, fall back to a subgradient” would make the wrong call on the second, which is why they are not one variant.

NonsmoothConePoint covers:

• the apex with a collapsed dual, and the boundary with a collapsed dual — the conic analogue of a weakly active row. Slack and multiplier vanish together, so the derivative is two-valued and depends on which way the perturbation pushes the block off its face. The NLP arm answers this class with a directional mode; the convex arm does not have one for cones yet, so it refuses rather than silently picking a side.
• a second-order boundary point too close to the apex, where w = (1, −s₁/s₀) would be built by dividing by round-off.
• a slack outside the cone beyond the solve’s own tolerance: there is no face to linearize against.
• a strictly interior block that does not complement (⟨s, z⟩ ≫ 0) — not a converged optimum, whatever its status field says.
• a PSD block where strict complementarity fails (rank Z ≠ n − rank S). That equality is what makes ker S the whole normal direction; without it a direction exists along which slack and multiplier vanish together, and dx/db is two-valued along it.
• the exponential and power cones’ degenerate faces (y = 0, z = 0), where the boundary has no tangent plane. There is deliberately no guard for the power cone’s |x| = 0 kink: x = 0 on the boundary forces y^α z^{1−α} = 0, i.e. one of those faces, so with y, z > 0 the two smooth sheets x = ±g never meet. A guard there would be unreachable code that reads like coverage.
• a non-symmetric dual off the facet’s normal ray. At a facet interior the normal cone is ℝ₊∇φ, so z = ν∇φ is the optimality condition, not an approximation.

ActiveSetOverdetermined has one case today:

• an apex-pinned block whose active set cannot absorb db. The apex is the one face that pins its whole block, so the step lives in ker(B) while feasibility needs A·dx = db. Where the two cannot both hold, what would come back is a least-squares compromise rather than a derivative. The model itself is usually perfectly smooth here: the guard fires where the classifier switched to Apex, and a decade further from the tip the boundary face returns the same derivative.

  The criterion is

  rank([A; B])  ==  rank(A) + rank(B)
  

  which is exact, not a dimension count. The quantity that matters is dim A(ker B) — the perturbations the step can actually reach — and the rank identity gives it as rank([A;B]) − rank(B); requiring that to be all of range(A) rearranges to the line above. An earlier dimension count n − rank(B) ≥ rank(A) is implied by it and strictly weaker: it passed a model whose equality lay entirely inside the pinned coordinates, where A(ker B) = {0} and no perturbation is reachable, and that model was served with an answer 33% off at every step size.

  rank(A), not A’s row count — a redundant equality does not shrink the space a step must reach. And B is the active rows that cannot be released: the cone faces and the active orthant rows. Active variable bounds are deliberately excluded, even though a bound pins its coordinate for the plain parametric_step — the release path can open a bound, and refusing at build time would take that path away too.

  One thing it is not: a promise that every db would have failed. A build serves every later perturbation and cannot know which are coming, so it refuses on the existence of one unreachable direction. That is deliberate, and stronger than “no answer exists here”.

  The complementary case is not a refusal at all. Ask a served build for a db outside range(A) — or take a plain step on a bound-pinned model the exclusion above deliberately serves — and the perturbed problem is simply infeasible: there is no derivative, and what comes back is a least-squares answer to an unanswerable question. ill_conditioned() is what tells you, after the step and never at build time. It is the residual clause that fires, never the condition estimate: the regularized KKT is perfectly well conditioned on these models (3.0e10 against a 1e14 threshold), so checking ill_conditioned() straight after build_conic returns false. Measured residuals are 0.333 and 0.8 against a 1e-6 threshold. Take the step, then check — and where a bound is what pins the model, parametric_step_bounded reproduces the re-solve exactly.

Two of these thresholds are calibrated against the non-symmetric driver, whose accuracy is well short of the symmetric IPM’s, and the measured populations are recorded at their definitions rather than left as round numbers — the first value tried for the dual-ray test refused two of four correct solutions.

The apex/boundary decision is relative to the problem’s primal scale (max(‖h‖∞, ‖Gx‖∞, 1)), the same quantity the orthant guard above uses, so the two cannot disagree about what “zero” means on one solution.

For the NLP arm’s much larger sensitivity surface — fix-relax and path modes, the directional decision at a kink, the corrector, activity classification, and the covariance/identifiability statistics — see Sensitivity Analysis. The two arms are not at parity today.

Presolve (PaPILO-inspired)

Before the interior-point solve, POUNCE can apply a transaction-stack presolve with full primal and dual postsolve, modeled on PaPILO. The catalog:

• empty / duplicate / parallel (scalar-multiple) rows,
• fixed-variable elimination (singleton equalities),
• free columns and free-column singletons,
• activity-based redundancy and infeasibility detection,
• forcing constraints (a row at its activity extreme pins its variables),
• dominated columns (sign-definite columns optimal at a bound),
• bound tightening (domain propagation), with the active-bound multiplier re-attributed to its source row in postsolve,

iterated to a fixpoint so reductions cascade. Each reduction carries the data to reverse itself, and the postsolve reconstructs a valid KKT point of the original problem — the dual recovery is the contract, and is verified by KKT-residual tests. A cone-aware variant (presolve_conic) gates the ≤-row reductions off second-order-cone blocks (which are coupled) and recovers the reduced cone partition.

The iteration also carries a layer cap, and on a model with a long bound-propagation chain — commonly, on roughly half the LP corpus — the cap is what stops it rather than the fixpoint. That distinction is visible: presolve reports which of the two happened and the CLI says so on its summary line (see LP / QP Solver Routing). A truncated reduction is still correct — every reduction it did apply is a sound transform with its own dual recovery — and measured across the LP and QP suites the truncation costs only box tightness, never a structural reduction.

Presolve is applied automatically on the CLI LP/QP route; it lives in pounce-convex::presolve for Rust callers. See LP / QP Solver Routing.

Differentiable convex layers (JAX)

pounce.jax exposes the solve as a differentiable JAX op via the implicit-function theorem on the KKT system at the optimum (Amos & Kolter, OptNet, 2017). The forward calls the solver; the backward is a single linear solve through the same KKT matrix.

import jax, jax.numpy as jnp
from pounce.jax import solve_qp, solve_socp, QpLayer

# x*(c) for a parametric QP, differentiable w.r.t. all of P, c, G, h, A, b.
def loss(c):
    x = solve_qp(P=P, c=c, G=G, h=h)
    return jnp.sum((x - target) ** 2)

grad_c = jax.grad(loss)(c0)        # exact gradient via implicit diff
J = jax.jacrev(lambda c: solve_qp(P=P, c=c, G=G, h=h))(c0)

• Gradients are provided w.r.t. every parameter that enters through the optimum: c, b, h, and the matrices P, G, A (the full OptNet matrix derivatives; ∇P is the symmetric gradient).
• solve_socp differentiates SOCPs too — the complementarity row uses the cones’ arrow operators in place of the orthant’s diagonal.
• QpLayer captures a fixed P/G/A structure for use inside a larger JAX model, with jax.grad / jacrev / vmap and a parallel .batch.
• A warm start may be passed through (non-differentiated — it cannot change the solution or its gradients, only the iteration count).

All gradients are validated against finite differences in the test suite.

Global Optimization

Most of POUNCE settles a problem at a local optimum (the NLP filter-IPM and SQP) or exploits convexity so that local is global (the convex/conic IPM). For a genuinely nonconvex problem, POUNCE offers one certified-global route, and it is for polynomials:

• The SOS / Lasserre hierarchy (pounce-convex) — for polynomial problems, via a single semidefinite program. Callable from Rust (sos_minimize) and Python (pounce.sos_minimize).

It returns a result that is certified: a lower bound together with a moment certificate that, when exact, pins the global minimum and recovers its minimizer(s).

There is no general-purpose spatial branch-and-bound solver in POUNCE. For a nonconvex problem that is not polynomial — anything with exp/ln/trig — POUNCE has no certified-global path. Use the local NLP solver from several starting points (see the multistart notebooks below), or reformulate into the convex cone library. A pounce-global crate was prototyped and removed from main before release; its design is recorded in dev-notes/spatial-bnb-design.md.

The SOS / Lasserre path (polynomials)

When the objective and constraints are polynomials, the sum-of-squares / moment approach in pounce-convex certifies the global minimum from a single semidefinite program — no branching — by searching for the largest γ such that p(x) − γ lies in the Putinar cone (a sum of squares plus constraint multipliers). The SDP is solved by POUNCE’s own convex conic interior-point method; flat truncation of the resulting moment matrix certifies when the bound is exact, and a facial-reduction step recovers every global minimizer — even when the optimum is attained at several points.

From Python, a polynomial is a dict mapping an exponent tuple to its coefficient (the all-zeros key is the constant term):

from pounce.sos import sos_minimize

# x**4 - 2 x**2 + 3  ->  global minimum 2, attained at BOTH x = +1 and x = -1
r = sos_minimize({(4,): 1.0, (2,): -2.0, (0,): 3.0})
r.lower_bound       # ≈ 2.0
r.is_exact          # True — flat-truncation certificate: the bound is the minimum
r.minimizers        # both x = +1 and x = -1

Constraints are polynomials too, passed as inequalities (g_i(x) ≥ 0) and equalities (h_j(x) = 0); raise the relaxation order to tighten the bound (the Lasserre hierarchy) at the cost of a larger SDP. A runnable walkthrough — double well, a constrained problem, and a 2-D example — is in 18_sos_global_optimization.ipynb.

The same solver from Rust, via the pounce-rs facade with the convex feature on (pounce-rs = { version = "0.9", features = ["convex"] }):

#![allow(unused)]
fn main() {
use pounce_rs::convex::{sos_minimize, PolyProblem, Polynomial};
use pounce_rs::linsol::backend;      // the sparse LDLᵀ factory the solver takes

// x⁴ − 2x² + 3 → global minimum 2 at x = ±1.
let p = Polynomial::new(1, vec![(vec![4], 1.0), (vec![2], -2.0), (vec![0], 3.0)]);
let sol = sos_minimize(&PolyProblem::new(p), None, backend);
// sol.lower_bound ≈ 2; when the moment matrix is flat, sol.minimizers holds
// the global minimizer(s) — here both x = +1 and x = −1.
}

The full treatment lives in the pounce_convex::sos module documentation — reachable without a second dependency, since pounce_rs::convex re-exports the pounce_convex crate itself for anything outside its curated surface.

When SOS fits: polynomials of modest degree and dimension — one SDP, recovers all global minimizers, but the SDP grows with the relaxation order.

When SOS does not fit

For a general factorable problem (exp/ln/trig), or a polynomial whose SDP would be too large, the textbook tool is spatial branch-and-bound — and POUNCE does not have one. Two things you can do:

• Reformulate into the cone library. If the model can be cast as an LP, convex QP, SOCP, or an exponential / power / PSD cone program, local is global and the guarantee comes for free. See Choosing a Solver.
• Multistart the local solver. Running the NLP filter-IPM from many starting points finds the low minima in practice, but certifies nothing — there is no bound proving you have the global one. Three notebooks work through the tactics: repulsion-based sampling (19_find_minima_repulsion.ipynb), random restarts (20_find_minima_restart.ipynb), and basin hopping (21_find_minima_hopping.ipynb).

Solution Output

The .sol file

Following the AMPL solver convention, solving a positional .nl file writes a sibling <stub>.sol next to it — pounce problem.nl produces problem.sol. The file carries the primal x and dual lambda blocks plus an objno line with the AMPL solve_result_num, so AMPL (or any .sol reader) can pull the solution back:

pounce problem.nl                       # writes problem.sol
pounce problem.nl --sol-output out.sol  # write to an explicit path
pounce problem.nl --no-sol              # skip the .sol write

A .sol is written even when the solve fails, so the solve_result_num is always recoverable. Built-in problems (--problem …) have no .nl stub, so they only produce a .sol when --sol-output is given explicitly.

Reading solve_result_num

The objno line carries an AMPL solve_result_num (Gay 2005, Hooking Your Solver to AMPL §5). Consumers key on the band, not the exact number:

Pyomo maps each band to a TerminationCondition, so anything in 200–299 arrives as TerminationCondition.infeasible.

Solved: strict, acceptable, and square

POUNCE writes the same codes IPOPT’s AMPL driver writes, so a model can be moved between the two solvers without a reader change:

All three are accepted solves and all three sit in the 0–99 band. That matters beyond tidiness: Pyomo’s legacy .sol reader loads the 0–99 band as SolverStatus.ok and the 100–199 band as SolverStatus.warning, with TerminationCondition.optimal either way. POUNCE reported acceptable-level solves as 100 up to and including 0.10.0, so Pyomo logged

WARNING: Loading a SolverResults object with a warning status into model...
    - termination condition: optimal
    - message from solver: POUNCE 0.10.0: SolvedToAcceptableLevel

on a solve IPOPT loads clean (#591). The distinction between strict and acceptable convergence is still there — in the code itself (1, not 0), in the .sol message line, and in the JSON report’s status field — it simply no longer reads as a warning.

Feasible_Point_Found was 100 up to and including 0.10.0 for the same reason and was fixed the same way (#815). It is worth being precise about why a feasible point counts as solved here, because in general it would not: POUNCE emits this status only when the problem is square, which is the condition IPOPT uses for its own 2. On a square problem there is nothing to optimise — the objective is constant over a feasible set the equalities have already pinned — so a point that satisfies the constraints is the solution, and there is no further criterion it could be said to have missed. For a non-square problem the status is never produced.

The stakes were higher than a logged warning. Pyomo’s newer .sol reader (pyomo.contrib.solver) maps the 100–199 band to TerminationCondition.error, not to optimal-with-a-warning, so on that route a square flowsheet that POUNCE had solved to a constraint violation of 2.2e-06 was delivered to the caller as a solver error.

Infeasible: proved vs. local

Within the infeasible band POUNCE distinguishes how it knows:

The difference is real, not cosmetic. 201 is a structural detection made on the model’s bounds before iterating, not a certified proof — it is subject to the same floating-point limits as any interval computation, and is withheld whenever the violation is smaller than the feasibility tolerance. 200 is different in kind — on a nonconvex problem a positive local minimum of the violation does not rule out a feasible point elsewhere, which is why the console message says “Problem may be infeasible.”

Because 200 is an inference rather than a proof, it is withdrawn when POUNCE holds a point that contradicts it. Before any numerical path reports 200, the model’s own starting point is evaluated against every constraint; if it satisfies them all, the feasible set is demonstrably non-empty and the verdict becomes Error_In_Step_Computation (500) — an honest “the solve broke down” rather than a wrong answer. This can only ever withdraw a verdict: a model with no feasible point cannot produce such a point, so a correct 200 is unaffected. Supplying a feasible starting point is therefore worth doing on a model you believe is feasible but POUNCE reports otherwise.

When the region is found empty the solve is skipped entirely and the message names how it was found, so the claim is checkable:

POUNCE 0.9.0: InfeasibleProblemDetected (detected by presolve: bound propagation)
objno 0 201

201 requires presolve to be enabled (presolve=yes); it is off by default. A presolve-derived infeasibility is only reported when the contradiction holds on the original box — one produced by presolve’s own auxiliary elimination is re-checked after rollback and never certified.

One more route to 200: over-determined systems

An over-determined model — more equality rows than free variables, such as x == 0.2 with x == 0.8 — cannot be solved at all: it fails a structural gate before the first iteration. That used to be reported as Not_Enough_Degrees_Of_Freedom (504, the failure band), which says “cannot attempt this” for a model whose answer is already decided.

POUNCE now checks such a model for a bound-propagation contradiction on that failure path and reports 200 when it finds one. This does not need presolve=yes — nothing is transformed and no solve runs through the check — so it is the one way to reach the infeasible band with the default options and no iterations. A consistent over-determined system is unaffected and still reports 504.

Because the solve provably cannot run here, this route measures constraint residuals against each row’s declared magnitude rather than an absolute tolerance, so the verdict does not change when every row is multiplied by a constant. Elsewhere — wherever a solve can run — an infeasibility smaller than the feasibility tolerance is still withheld, as described above.

Choosing an output format

The .sol and JSON outputs are not exclusive — you can request both in the same run.

JSON Solve Report

Pass --json-output PATH to write a structured solve report alongside the regular console output:

pounce problem.nl --json-output result.json
pounce problem.nl --json-output result.json --json-detail full

The report carries everything an AMPL .sol file holds — status, primal x, dual lambda, suffix blocks — plus FAIR-aligned provenance metadata (Wilkinson et al. 2016, DOI 10.1038/sdata.2016.18) and, optionally, the per-iteration trajectory.

Detail levels

Choose summary for production logs and batch runs; full for debugging — it is the JSON equivalent of upstream’s print_level=8.

Reproducibility: recorded environment overrides

Solve-affecting environment variables — the POUNCE_FERAL_* linear-solver knobs and the legacy FERAL_PIVTOL / FERAL_PARALLEL — are captured into fair_metadata.environment when set, so a run that differs because one was exported in a shell profile says so instead of differing silently:

"environment": [
  { "name": "POUNCE_FERAL_PIVTOL", "value": "1e-6" }
]

The block is omitted entirely when no such variable is set (the common case). Debug-only gates (POUNCE_DBG_*) are not captured. See the schema reference for the full field contract.

Schema stability

The schema is versioned (pounce.solve-report/v1) so downstream tooling can pin against a major version:

• Adding fields is non-breaking — consumers must tolerate unknown fields.
• Removing or renaming a field bumps the major version (v1 → v2).

The Schema v1 Reference documents every field, the FAIR mapping, and the versioning policy in full.

POUNCE solve-report schema, v1

Schema tag: pounce.solve-report/v1

This document is the canonical reference for the JSON solve report emitted by pounce --json-output PATH and pounce_sens --json-output PATH. The report carries everything an AMPL .sol file holds — status, primal x, dual lambda, suffix blocks — plus FAIR-aligned provenance metadata and (optionally) the per-iteration trajectory.

Implementation: the serde structs live in crates/pounce-solve-report/src/lib.rs (per-iteration IterRecord in crates/pounce-nlp/src/solve_statistics.rs); crates/pounce-cli/src/solve_report.rs wires them to the CLI.

Why a structured solve report?

Production NLP workflows often need to (a) capture which solve produced which numbers for audit / reproducibility, (b) feed solver output into downstream tooling (notebooks, dashboards, ML pipelines) that don’t want to parse a free-form .sol file, and (c) compare runs across versions of pounce. Both upstream Ipopt’s stdout summary and AMPL’s .sol were designed for human consumption and AMPL’s reader respectively — neither carries provenance metadata, neither is schema-versioned, and neither is trivially machine-parseable across ecosystems.

A versioned JSON schema with FAIR-aligned provenance solves all three.

FAIR alignment

The fair_metadata block maps onto the four FAIR principles (Wilkinson et al. 2016, “The FAIR Guiding Principles for scientific data management and stewardship”, Scientific Data 3, 160018, DOI 10.1038/sdata.2016.18; citation verified via Crossref on 2026-05-14):

Versioning policy

schema is the version tag. Compatibility rules:

• Adding fields is non-breaking. Consumers MUST tolerate unknown fields. New optional fields land between versions; the major version doesn’t bump.
• Removing or renaming fields bumps the major version (v1 → v2). Consumers should pin against a major version (schema starts_with "pounce.solve-report/v1").
• Changing field semantics without a rename is forbidden. If semantics need to change, add a new field and deprecate the old.

The pre-1.0 phase of POUNCE itself does NOT relax this rule for the schema. Once a solve-report version ships, its field set is frozen even while the rest of the solver is under churn.

Top-level shape

{
  "schema": "pounce.solve-report/v1",
  "fair_metadata": { ... },
  "problem":       { ... },
  "solution":      { ... },
  "statistics":    { ... },
  "iterations":    [ ... ],  // optional, omitted when empty
  "linear_solver": { ... }   // optional, omitted when backend did not report
}

Fields

schema (string, required)

Identifier for this schema version. Always "pounce.solve-report/v1" for v1. Major-version bumps change the prefix; minor / patch (additive) changes do not.

fair_metadata (object, required)

solver sub-object

Input descriptor (input)

Tagged enum keyed on kind. Possible shapes:

{ "kind": "nl-file", "path": "/path/to/foo.nl", "size_bytes": 366 }
{ "kind": "builtin", "name": "rosenbrock" }
{ "kind": "tnlp-direct" }

• nl-file — the input came from .nl file at path. size_bytes is present when the file’s metadata is readable; consumers that want bit-exact provenance can hash the file themselves.
• builtin — the input was a built-in problem named by name (e.g. pounce --problem rosenbrock).
• tnlp-direct — used by library callers building a TNLP in-process without a .nl round-trip.

Environment overrides (environment)

An array of { "name", "value" } objects, one per solve-affecting environment variable set in the process at report time:

"environment": [
  { "name": "POUNCE_FERAL_PIVTOL", "value": "1e-6" }
]

The whole array is omitted when none are set (the common case). Only the variables that change pounce’s numerics or parallelism are captured — the POUNCE_FERAL_* linear-solver knobs and the legacy FERAL_PIVTOL / FERAL_PARALLEL. These alter the factorization and can otherwise silently differ a run between two machines (e.g. one with POUNCE_FERAL_PIVTOL exported in a shell profile) with nothing in the report saying so. The POUNCE_DBG_* debug gates are deliberately not captured — they only add diagnostic output and never change the result.

Presence records that the variable was set, not that it took effect: an explicit OptionsList setting (e.g. feral_pivtol in an options file) takes precedence over the env fallback. See Options › Environment overrides for the option each variable maps to.

problem (object, required)

Problem dimensions reported by the TNLP at get_nlp_info().

solution (object, required)

Suffix entries

{
  "name": "sens_sol_state_1",
  "target": "var",
  "kind": "real",
  "values": [0.576..., 0.378..., -0.046..., 4.5, 1.0]
}

statistics (object, required)

Projection of pounce_nlp::solve_statistics::SolveStatistics minus the per-iteration history (which lives at the top level when present).

null values. The four residuals are produced by the convergence check at the end of a solve. A solve the solver refused — rejected during setup (NotEnoughDegreesOfFreedom, InvalidProblemDefinition), aborted, or caught by the batch panic handler — never reaches it, and these slots are emitted as null rather than 0.0. A zero there is indistinguishable from a perfect solve, and consumers acted on it: it was enough to make pounce.minimize report success=True for a problem the solver had declined to attempt.

The two objective fields follow the same rule for the same reason: 0.0 is an ordinary objective value, so it cannot signal “never evaluated”. They are seeded from the current iterate whenever one exists, so they are null only when the solve produced no point at all.

Consumers should treat null as “not computed”, not as zero. pounce’s own readers map it to NaN, which fails closed against any value <= tol test. | num_obj_evals | integer | eval_f call count. | | num_constr_evals | integer | eval_g call count. | | num_obj_grad_evals | integer | eval_grad_f count. | | num_constr_jac_evals | integer | eval_jac_g count. | | num_hess_evals | integer | eval_h count. | | total_wallclock_time_secs | float | Wall time spent inside optimize_*. | | restoration_calls | integer | Number of restoration-phase entries (pounce#12). | | restoration_inner_iters | integer | Cumulative inner-IPM iterations across all restoration calls. | | restoration_outer_iters | integer | Outer iterations that ran in restoration mode (R-line equivalents). | | restoration_wall_secs | float | Wall time spent inside perform_restoration. | | quality_escalations | integer | Times the linear solver escalated its factorization (IncreaseQuality) during this solve, restoration sub-solves included (pounce#857). 0 on every backend that cannot escalate and on any solve that never stalled. Present since the field was added; reads 0 on reports written before it, since it deserializes with a default. It is the only trace an escalation leaves — status, objective, iteration count and engine are all unchanged by one, and the console’s q info-string flag misses every escalation taken inside restoration, because those rows carry no info column. On a laddered run this is the promoted solve’s count, not the base solve’s — the same rule iteration_count follows. It is a sharp edge here because feral_increase_quality_retry promotes a re-solve that by construction escalated zero times, so a run whose base solve escalated twenty-five times reports 0 once the recovery lands; second_opinion.base_status records what it recovered from, and feral_increase_quality_retry=no reproduces the base solve outright. |

Eval counters (num_*_evals) populate only on the .nl-file path because the pounce binary’s CountingTnlp wrapper tracks them. Library callers using IpoptApplication::optimize_tnlp directly see zeros there; the underlying counts are still available through upstream’s IpoptCalculatedQuantities if needed.

iterations (array of object, optional)

Per-iteration trajectory. Emitted only at --json-detail full (when IpoptApplication::enable_iter_history() was called). Omitted from JSON entirely when empty.

Each row maps to one line of the upstream-formatted console iter table. Fields:

linear_solver (object, optional)

Aggregate post-mortem from the symmetric-indefinite linear backend that solved the KKT systems. Populated only when the backend self-instruments (the default FERAL backend does; HSL MA57 and custom backends plugged through set_linear_backend_factory do not). Omitted from JSON when no backend reported.

Detail levels

The --json-detail LEVEL flag selects how much detail is emitted. Levels map to verbosity in the same spirit as upstream’s print_level (0 silent → 12 maximum debug):

summary is the right choice for production logs and batch runs. full is the debugging equivalent of upstream’s print_level=8.

Worked example

pounce_sens crates/pounce-cli/tests/fixtures/parametric.nl out.sol --json-output result.json --json-detail full produces (truncated for brevity):

{
  "schema": "pounce.solve-report/v1",
  "fair_metadata": {
    "result_id": "1778777029606881000-76543",
    "created_at_iso": "2026-05-14T16:43:49.606Z",
    "created_at_unix_nanos": 1778777029606881000,
    "elapsed_seconds": 0.011,
    "solver": {
      "name": "pounce",
      "version": "0.1.0",
      "target_triple": "x86_64-apple-darwin"
    },
    "license": "EPL-2.0",
    "input": {
      "kind": "nl-file",
      "path": "crates/pounce-cli/tests/fixtures/parametric.nl",
      "size_bytes": 366
    }
  },
  "problem": { "n_variables": 5, "n_constraints": 4, "n_objectives": 1, "minimize": true },
  "solution": {
    "status": "SolveSucceeded",
    "status_upstream": "Solve_Succeeded",
    "solve_result_num": 0,
    "objective": 0.5510204081632656,
    "x":      [0.6326530575201161, 0.3877551079678144, 0.020408165487930466, 5.0, 1.0],
    "lambda": [-0.16326530000405073, -0.28571431357898697, -0.16326530000405073, 0.18075803406303625],
    "suffixes": [{
      "name": "sens_sol_state_1",
      "target": "var",
      "kind": "real",
      "values": [0.5765305974643309, 0.3775510440570709, -0.04591835847859835, 4.5, 1.0]
    }]
  },
  "statistics": { "iteration_count": 9, "final_dual_inf": 2.89e-14, "...": "..." },
  "iterations": [
    { "iter": 0, "objective": 0.0451, "inf_pr": 5.0, "inf_du": 0.407, "mu": 0.1,
      "d_norm": 0.0, "regularization": 0.0, "alpha_dual": 0.0, "alpha_primal": 0.0,
      "alpha_primal_char": " ", "ls_trials": 0 },
    { "iter": 1, "objective": 0.957, "inf_pr": 0.212, "...": "..." }
  ]
}

Consumer guidance

• Pin the major version. Check schema.startswith("pounce.solve-report/v1") before consuming.
• Tolerate unknown fields. New optional fields will land between minor versions of pounce. Use serde(default) / equivalent.
• Distinguish “no solve” from “solve produced zero”. Pre-solve, scalar fields are 0.0 (not NaN, because JSON has no NaN literal). statistics.iteration_count == 0 is the signal that no solve occurred.
• solution.x / solution.lambda may be empty. When the binary couldn’t capture the final iterate (currently: the pounce binary on its newton_driver fast-path for m=0, n≤1000 problems), the arrays are empty and the keys are omitted from JSON entirely. pounce_sens always populates them.

References

• Wilkinson et al. (2016). “The FAIR Guiding Principles for scientific data management and stewardship.” Scientific Data 3, 160018. DOI 10.1038/sdata.2016.18. (Verified via Crossref 2026-05-14.)
• Gay (2005). “Hooking Your Solver to AMPL.” https://ampl.com/REFS/hooking2.pdf. §5 (Returning Results to AMPL) for the .sol baseline this schema is structured around.
• SPDX license identifiers: https://spdx.org/licenses/.

Verifying Solutions

pounce verify <problem.nl> <claim.sol> [OPTIONS]

pounce verify independently checks that the solution in a .sol file actually satisfies the constraints and bounds of a .nl problem. It re-derives feasibility from the model itself — it does not trust the .sol’s status line, and it does not rerun the solver. This makes it the trust anchor when pounce is a tool an agent calls: the agent proposes a solution, and a small, deterministic checker disposes.

Optimization is unusually well-suited to this because a solution is far cheaper to verify than to produce: a claimed x* is just numbers, and feasibility is a single constraint evaluation — g_l ≤ g(x*) ≤ g_u, x_l ≤ x* ≤ x_u — O(nnz) work with no resolve and no dense linear algebra.

Status. The verify check itself — recompute feasibility against the canonical model, with a content-addressed receipt — is solid and ready to use; it needs no secrets and is the recommended default. The signing and remote-service trust layer layered on top of it (HMAC receipts, the signer_service.py reference, running the MCP server as a remote authority) is a proof of concept: it demonstrates the architecture but is not hardened for production. If you want to rely on the signed/remote path for real, see Status and hardening at the end for the checklist of what that would take.

What it defends against

In an agent workflow, three things can go wrong with “here is a solution”:

The key design rule: always verify against the canonical problem, never against whatever the agent claims it solved. If the agent loosened a bound to manufacture feasibility, the returned x* still violates the canonical bound, and verify reports it.

Output and exit codes

$ pounce verify gaslib40_steady.nl good.sol
pounce verify — independent solution check
  problem : gaslib40_steady.nl  (1694 vars, 1682 cons)
            sha256:4bb435a3…
  solution: good.sol
            sha256:b77d9e7b…
  claimed solve_result_num: 0

  feasibility (tol 1.0e-6):
    max constraint violation: 1.407e-12  at c[114] (value 1.4e-12, bounds [0, 0])
    max bound violation     : 9.775e-9   at x[24]  (value 1.05, bounds [1.05, 2.0])
  objective at x*: 1.2899875310e0

  optimality (tol 1.0e-6, duals + bound multipliers supplied):
    KKT stationarity residual (bound-projected)  : 2.675e-3  (dual sign +1)
    dual infeasibility (with z_L/z_U suffixes)   : 6.248e-14
    constraint complementarity (rows, |λ|·slack) : 0.000e0
    bound complementarity (vars, |z|·slack)      : 9.091e-10

  VERDICT: VERIFIED — solution is feasible for the canonical problem

A consumer (CI step, agent harness, Makefile) gates on the exit code.

Options

Feasibility is the gate; optimality is reported

By default only feasibility gates the exit code. Feasibility is rigorous and sign-convention-independent — it is the guarantee that matters when the claim is “this solution meets the constraints.”

When the .sol carries constraint duals, verify also reports a KKT stationarity residual (the bound-projected “dual infeasibility”: the part of ∇f + Jᵀλ that a valid sign-constrained bound multiplier cannot absorb) and the complementarity residuals below. These are informational unless you pass --require-optimal. The AMPL dual-sign convention can differ from pounce’s, so verify computes the residual for both signs and reports the better one plus the sign it used.

The two complementarity residuals

Two distinct quantities answer to “complementarity”, and at the same point they can disagree by many orders of magnitude. verify names both by what they range over, and never prints an unqualified complementarity residual:

The bound one is what a solver prints as Complementarity — Ipopt’s and pounce’s own end-of-solve report alike. Do not compare a solver’s Complementarity against the row line; they measure different things and neither is wrong for what it measures.

When the .sol carries no ipopt_zL_out / ipopt_zU_out suffixes, bound complementarity is reported as not checked, never as 0.0.

Bound multipliers sharpen the stationarity check

Without the suffixes, z_L and z_U are inferred from which bounds are active, and the reported stationarity residual is bound-projected: it projects out exactly the component a bound multiplier would carry, so it reads 0.0 on a point whose bound multiplier is missing or wrong.

When the .sol does carry them (pounce always writes them; so does Ipopt’s AMPL interface), verify additionally reports the exact dual infeasibility ‖∇f + Jᵀλ − (zL_out + zU_out)‖∞, using the multipliers the file actually claims. That number is directly comparable to the solver’s Dual infeasibility, and it is what --require-optimal gates on whenever it is available — the projected residual can only understate it.

The JSON receipt

--json-output writes a machine-readable receipt that content-addresses both inputs by SHA-256 — so a downstream consumer can confirm exactly which problem and which solution were checked:

{
  "pounce_verify_version": 1,
  "solver": "pounce 0.4.0",
  "problem":  { "path": "…", "sha256": "4bb435a3…", "n_vars": 1694, "n_cons": 1682 },
  "solution": { "path": "…", "sha256": "b77d9e7b…", "duals_present": true },
  "tolerances": { "feasibility": 1e-6, "optimality": 1e-6 },
  "feasibility": {
    "max_constraint_violation": 1.4e-12,
    "worst_constraint": { "index": 114, "name": "c[114]", "value": 1.4e-12,
                          "lower": 0.0, "upper": 0.0, "violation": 1.4e-12 },
    "max_bound_violation": 9.77e-9,
    "worst_bound": { "index": 24, "name": "x[24]", … },
    "feasible": true
  },
  "optimality": {
    "available": true,
    "stationarity_residual": 2.6e-3,
    "stationarity_residual_with_bound_multipliers": 6.2e-14,
    "constraint_complementarity_residual": 0.0,
    "bound_complementarity_residual": 9.1e-10,
    "bound_multipliers_present": true,
    "complementarity_residual": 0.0,
    "optimal": true, "note": "…"
  },
  "verdict": "VERIFIED",
  "verified": true
}

bound_complementarity_residual and stationarity_residual_with_bound_multipliers are null when the .sol carries no ipopt_zL_out / ipopt_zU_out suffixes — null means not checked, not zero. complementarity_residual is a deprecated alias of constraint_complementarity_residual, kept so v1 consumers keep parsing; its bare name is the one that invited the wrong comparison, so read the qualified field instead.

A consumer should accept a solution iff:

1. verified == true, and
2. problem.sha256 equals the SHA-256 of its own canonical .nl, and
3. (when signing is used) the signature validates — see below.

Checking the hash in step 2 is what closes the “solved the wrong problem” gap at the receipt layer: the receipt is only meaningful for the exact problem bytes it names.

The default: recompute, don’t trust a receipt

The strongest and simplest design uses no key and no signature at all: the consumer runs pounce verify itself, against its own copy of the canonical .nl.

# the consumer does this — not the agent
pounce verify ./canonical/problem.nl ./from-agent/claim.sol || reject

Because verification is keyless, deterministic, and cheap (O(nnz), no resolve), the consumer can afford to just do it rather than trust someone else’s word. In this design the agent is never in the trust path: it hands over x*, and the consumer believes its own arithmetic. There is no key to steal, so the question “what if the agent gets the key?” does not arise. Forgery is impossible because nothing is being trusted on faith — feasibility is decided by evaluating g(x*), not by matching fields in a document.

This is the recommended default. Prefer it whenever the consumer can run pounce verify (or call a verifier it controls). Reach for signatures only when it genuinely cannot — see below.

Signed receipts — trust transport, conditional on key isolation

Signing addresses a narrower situation: the consumer won’t or can’t recompute — a remote or expensive verifier, or an audit log you want to trust later without re-solving — and instead wants to trust a receipt produced elsewhere. A signature lets that receipt be checked without redoing the work.

When the POUNCE_VERIFY_KEY environment variable is set (non-empty), the receipt gains:

"signature_alg": "HMAC-SHA256",
"signed_fields": ["verify_version","nl_sha256","sol_sha256",
                  "n_vars","n_cons","feasible","verified","verdict"],
"signature": "5bdcc146bf60754e…"

The signature is HMAC-SHA256(key, preimage), where preimage is a deliberately float-free byte string — only hex hashes, integer counts, and the verdict — so any language reproduces it byte-for-byte without float-formatting parity problems. The exact preimage is:

pounce-verify-receipt/v1
verify_version=1
nl_sha256=<hex>
sol_sha256=<hex>
n_vars=<int>
n_cons=<int>
feasible=<true|false>
verified=<true|false>
verdict=<VERIFIED|REJECTED>

(eight lines, \n-joined, with a trailing newline; booleans lowercase.) A holder of the key recomputes the HMAC over this preimage and compares it to signature.

What the signature does and does not guarantee

HMAC gives existential unforgeability under chosen-message attack — but only while the key stays secret. That single condition carries the entire weight of the scheme:

• A party without the key cannot mint a receipt that validates. It cannot produce a correct MAC for any message it has not already seen one for. This is the real guarantee.
• A party with the key can sign anything — including a hand-written receipt that falsely says verified: true. If the agent obtains the key, signing buys you nothing, and you must fall back to recomputation.

So the security reduces to a question that is not cryptographic: can the agent reach the key? If the keyholder and the agent share a process, a user account, an environment, or a filesystem, the agent simply reads it — echo "$POUNCE_VERIFY_KEY", cat keyfile, or introspecting shared memory — and the signature becomes theater. An agent with a shell on the same user/host as the signer is not kept out by HMAC.

Setting POUNCE_VERIFY_KEY for a pounce verify that runs in the agent’s own environment is the weakest form and provides no real protection against that agent. It is only meaningful when the key lives behind a privilege boundary the agent cannot cross. See Out-of-process signing below.

Two further limits HMAC does not address, even with a perfectly isolated key:

• Replay. A previously, legitimately-signed receipt can be re-presented. Here that is benign: a receipt attests “this specific x* (by sol_sha256) is feasible for this specific problem (by nl_sha256),” and that statement stays true. Replay cannot manufacture a valid-looking attestation for a new (problem, solution) pair the signer never blessed. If you ever need freshness (“the work was done this time”), add a nonce/timestamp to the signed preimage — but for “does this solution meet the constraints,” freshness is irrelevant.
• Wrong model. The signature says nothing about whether the .nl encodes the problem you meant. That spec-ownership question must live outside the agent regardless of any crypto. verify certifies feasibility against a given model; it does not certify that the model is the right one. Likewise, feasibility is fully checkable but global optimality is not — the stationarity residual certifies a first-order/KKT point, not a global minimum.

Out-of-process signing

For a signature to mean anything against a capable agent, the key must live in a separate trust boundary — a different process and user/container/ host that the agent cannot exec into or read. The signer exposes “verify-and-sign” as an RPC; the key never enters the agent’s address space.

 ┌─ agent boundary (no key) ────────────┐      ┌─ signer boundary (holds key) ─────┐
 │ agent proposes x*                     │      │ POUNCE_VERIFY_KEY in THIS env only │
 │   │                                   │ RPC  │   1. run `pounce verify` on the    │
 │   └── POST /verify {nl, sol} ─────────┼─────▶│      CANONICAL .nl (+ the key)     │
 │                                       │      │   2. binary signs the receipt      │
 │   signed receipt ◀────────────────────┼──────┤   3. return receipt JSON           │
 └───────────────────────────────────────┘      └────────────────────────────────────┘
        │
        └── relays receipt to the consumer
consumer: accept iff  verified==true  ∧  problem.sha256==canonical  ∧  signature valid

What each party can do under this split:

The boundary is only real if the agent cannot run code as the signer’s user or on its host. Running the signer as a separate user, container, or host (or behind a KMS/HSM that signs without exposing the key) is what turns “signed” from theater into a guarantee. An MCP server is already a separate process from the model, which helps — but only achieves isolation if the agent also lacks a shell on the same user/host.

A minimal reference signer is in studio/mcp/examples/signer_service.py: a stdlib HTTP service that holds the key in its own environment, shells out to pounce verify, and returns the signed receipt. The agent calls it; the agent’s environment never contains the key.

Use in an agent workflow

Putting it together — recompute by default, sign only to transport trust:

agent ── proposes x* ──▶  consumer / verifier-it-controls
                            1. pin + hash the canonical .nl
                            2. pounce verify .nl .sol   (against the CANONICAL model)
                          ◀─ accept iff verified==true ∧ problem.sha256==canonical

When the verifier must be remote and the consumer won’t recompute, insert an out-of-process signer (above) and add ∧ signature valid to the consumer’s acceptance test — remembering that the last clause is only as strong as the signer’s key isolation.

The pounce-studio MCP server exposes verify_solution so an agent can request a check but cannot fake its result. Deploy that server as a distinct boundary from the agent (separate user/container) for the signature to carry weight; otherwise rely on the consumer recomputing.

Status and hardening

What is ready to use as-is:

• The feasibility check (pounce verify, and the consumer-recomputes pattern). It is deterministic, keyless, content-addressed, and rigorous — this is the part to build on.

What is a proof of concept — demonstrates the shape, not hardened:

• HMAC signing via POUNCE_VERIFY_KEY, the signer_service.py reference, and treating a remotely-deployed MCP server as a signing authority.

If you ever want to depend on the signed/remote path in production, these are the gaps to close. None are implemented here.

Key management

• Don’t keep the key in a plain environment variable or file. Use a KMS/HSM (or sealed secret) that signs without exposing the key to the process — then even a compromised signer host can’t exfiltrate it.
• Add key rotation and a key id in the receipt (kid) so a consumer knows which key to check against and old receipts stay verifiable across rotations.
• Consider an asymmetric scheme (e.g. Ed25519) instead of HMAC when more than one party must verify without also being able to sign — HMAC’s symmetric key means every verifier is also a forger. Public-key signatures give public verifiability with a single private signer.

Transport / service (the moment it leaves stdio)

• TLS on the endpoint; never plaintext for a service that holds a key.
• Authn/authz — bearer token or OAuth on every request (MCP’s HTTP transport supports this). An unauthenticated endpoint that runs solves and shells out is effectively remote code execution.
• Resource limits — request-size caps, solve timeouts (there is a timeout_seconds, but also wall/CPU/memory limits at the OS level), concurrency caps, and rate limiting.
• Sandbox the solve — treat every .nl as untrusted input. Parsing and evaluating an arbitrary model is attacker-controlled computation; run it in a locked-down container/user with no network and a constrained filesystem.

Input handling

• Over a network the path-based tools (nl_file/sol_file) assume a shared filesystem. Prefer content upload (the server receives and hashes the exact .nl/.sol bytes) so there’s no path-traversal surface and the receipt binds what was actually sent. The reference signer’s POUNCE_SIGNER_ROOT allowlist is a stopgap, not a substitute.

Freshness / replay

• The current preimage has no nonce or timestamp, so a signed receipt is replayable. That is benign for “is this x* feasible” (a timeless fact), but if a consumer needs “this was checked recently” or “in response to my request,” add a nonce/timestamp (and a receipt expiry) to the signed preimage and bump the pounce-verify-receipt version.

Auditability

• Log every verification (problem hash, solution hash, verdict, key id, caller identity) to an append-only store, so a disputed result can be reconstructed. Keep the key out of the logs.

Standing non-goals (true regardless of hardening)

• verify certifies feasibility against a given model — it does not certify the model is the right one. Model/spec correctness must be owned outside the agent.
• Feasibility is fully checkable; global optimality is not. The stationarity residual certifies a first-order/KKT point, not a global minimum.

Sensitivity Analysis

POUNCE includes a parametric sensitivity capability compatible with upstream Ipopt’s contrib/sIPOPT/ (Pirnay, López-Negrete & Biegler 2012, DOI 10.1007/s12532-012-0043-2). It computes the first-order change in the optimal primal solution with respect to a problem parameter, reusing the KKT factorization from the converged solve. Four entry points cover the common workflows.

AMPL CLI

The main pounce driver auto-detects the sIPOPT suffixes (sens_state_1, sens_state_value_1, sens_init_constr) in an input .nl, runs a post-optimal sensitivity step after the solve, and writes the perturbed primal back as a sens_sol_state_1 suffix — no separate binary or flag needed:

pounce problem.nl                   # writes problem.sol
pounce problem.nl out.sol --json-output result.json --json-detail full

pounce_sens is retained as a thin backward-compatibility alias: pounce_sens in.nl out.sol is identical to pounce in.nl out.sol, so existing AMPL / solver scripts keep working unchanged.

Related flags:

• --sens-boundcheck / --sens-bound-eps EPS — hold the perturbed primal x* + Δx at the declared bounds by pinning each crossing coordinate there and re-solving, so the others move with it.
• --compute-red-hessian / --rh-eigendecomp — compute the reduced Hessian (and its eigendecomposition) over the variables tagged by the red_hessian integer var-suffix.

The sIPOPT option names

The same requests can be made by upstream sIPOPT’s own option names — on the command line as key=value, or in an ipopt.opt — so a script written for sIPOPT keeps working. They are read by the pounce driver and by the pounce-sensitivity builder alike:

Two of these deliberately differ from upstream’s registered default, because honouring that default would change results for anyone who never set the option. run_sens is registered no upstream, but pounce runs the step whenever the input declares the suffixes, and only an explicit run_sens=no turns it off. sens_max_pdpert is registered 1e-3, but pounce applies no cap unless you set one — an unset sens_max_pdpert reports the step however hard the factor was regularized, as it always has. Check SensResult::kkt_perturbations (Python: info["kkt_perturbations"]) if you want to see the perturbation without capping on it.

n_sens_steps is the one sIPOPT key pounce does not honour: only the single sens_state_1 perturbation tier is implemented, so any value other than the default 1 is refused with an explanation rather than silently rounded down (gh#677).

Rust library

Reach the sensitivity path through the pounce-rs facade, with the sensitivity feature on:

[dependencies]
pounce-rs = { version = "0.9", features = ["sensitivity"] }

SensSolve is a builder that wraps the on_converged callback plumbing into a single call:

#![allow(unused)]
fn main() {
use pounce_rs::sensitivity::SensSolve;

let result = SensSolve::new(vec![2, 3])
    .with_deltas(vec![0.05, 0.0])
    .with_reduced_hessian()
    .run(&mut app, tnlp);
// result.dx, result.reduced_hessian, result.status
}

with_reduced_hessian_eigen() adds the eigendecomposition, and with_boundcheck(eps) enables the bound refinement described under Bending the estimate around a bound.

Eigenvector sign convention

Every eigendecomposition POUNCE hands back — the reduced Hessian’s here and through the CLI and Python wrappers, the QP one from QpSensitivity.reduced_hessian, and sens_covariance().eigen() / sens_information().eigen() in pyomo-pounce — returns sign-pinned eigenvectors: the largest-magnitude component of each column is positive, ties broken by the earliest row. v and -v are equally valid eigenvectors, so without a convention the direction you read back depends on the arithmetic that produced it and is not reproducible across builds or machines.

The sign is all that is pinned. A repeated eigenvalue leaves the basis within its eigenspace arbitrary — any rotation of those columns diagonalizes equally well — so read a degenerate block as a subspace, not column by column.

Python

solve_with_sens exposes the same capability from the cyipopt-compatible Python wrapper:

# pin_constraint_indices is required; pass deltas=..., compute_reduced_hessian=True,
# or both. Returns (x, info) — sensitivity outputs live in the info dict.
x, info = prob.solve_with_sens(x0, pin_constraint_indices=[2, 3],
                               deltas=[0.05, 0.0], sens_boundcheck=True)
# info["dx"], info["reduced_hessian"], info["reduced_hessian_eigenvalues"], ...

compute_reduced_hessian=True returns the reduced Hessian in info["reduced_hessian"]; rh_eigendecomp=True adds its eigendecomposition; sens_bound_eps=… tunes the bound refinement. See python/notebooks/04_sensitivity.ipynb for a walkthrough.

The analysis layer: pounce.sensitivity

solve_with_sens answers one perturbation. The pounce.sensitivity subpackage is the whole analysis surface — the step in every mode, what the step did about the bounds, the active-set events along a path, and the estimation statistics — over a session that is a solved NL plus the factorization it left behind:

import pounce
from pounce.sensitivity import solve_for_sensitivity, solution, covariance

nl = pounce.read_nl("model.nl")
sess = solve_for_sensitivity(nl, pins={"p": 4})     # 4 is p's pin row

solution(sess, [4], [0.05])                         # the moved solution
solution(sess, [4], [0.05], mode="fix_relax")       # bounds respected
solution_report(sess, [4], [0.05])                  # what it did about them
active_set_changes(sess, [4], [0.05])               # the events, in order

Estimation statistics need the fitted parameters’ .col columns and the residual variables’ columns, which is what declare_sens_fitted and declare_sens_residual resolve to on the Pyomo side:

sess = solve_for_sensitivity(nl, fit_rows={"a": 0}, res_rows={None: [1, 2, 3]})
cov = covariance(sess)          # cov["a"], cov.std_err["a"], cov.eigen()
inf = information(sess)         # the reduced Hessian it inverts

A parameter is addressed by the full-g row of the defining equality it is pinned by; a fitted parameter or residual by its full-x column. Keys are opaque — they order the result and label it — so results come back keyed by whatever the caller keyed the session with (names here, Pyomo component data under pyomo_pounce).

Two index-space warnings that a caller working in raw rows owns, and that a modelling layer would otherwise hide:

• .col order is full-x; the factor’s x block is var-x, which drops every variable the solve removed as fixed (lb == ub under the default fixed_variable_treatment=make_parameter). Route every factor index through session.primal_row(), which raises rather than returning a neighbouring variable’s answer (gh#450).
• fit_rows and res_rows hold variable columns, not constraint rows, including for the residuals.

pyomo_pounce is a caller of this package, not a reimplementation of it: sens_solution, sens_solution_report, sens_active_set_changes, sens_covariance and sens_information resolve Pyomo components to rows and hand the same session to the same functions. A fix here reaches both.

Pyomo

pyomo_pounce wraps the same machinery in a declare-then-query interface: flag the parameters that matter while building the model (no perturbed values required), solve normally, then ask for derivatives. Parameters are declared with declare_sens_param (mutable Param or fixed Var, scalar or indexed); when declarations are present, SolverFactory("pounce").solve(m) runs in-process and keeps the converged KKT factorization, so every query afterwards is a single backsolve.

A declared Param should enter the model through one defining equality: a single variable equal to the param, the shape a parameterized initial condition already has (m.x0 == m.x0_hat). Such a model solves as written, on every solve, and the defining equality is the row the machinery perturbs. A declared Param without that form is rewritten in place once, at declaration, with a warning: its occurrences are replaced by a substituted variable held by a new defining equality, the affected rows edited in place so their names are untouched. Writing the defining equality yourself avoids the rewrite and is the recommended form. Editing the model after declaration so a declared Param appears in new expressions is unsupported: re-declare on the current model instead. Repeated solves of one declared model, the receding-horizon pattern, pay no per-solve model copy and no per-solve rewrite.

import pyomo.environ as pyo
import pyomo_pounce
from pyomo_pounce import declare_sens_param, sens_jacobian, sens_solution

m.p = pyo.Param(initialize=2.0, mutable=True)
declare_sens_param(m.p)                 # a flag, not a perturbation

pyo.SolverFactory("pounce").solve(m)    # ordinary solve

sens_jacobian(m.x, wrt=m.p)                  # dx*/dp (float)
sens_jacobian(m.con, wrt=m.p)                # d(multiplier of con)/dp
sens_jacobian(m.obj, wrt=m.p)                # df/dp, the total derivative
G = sens_jacobian(m.z, wrt=m.r)              # containers -> Jacobian object
G[m.z[1], m.r[2]]; G.to_dataframe()     # element access / full Jacobian
sens_solution(m, [(m.p, 2.5)])               # first-order solution estimate at
                                        # new values, clamped to bounds

The objective: df/dp

of= the model’s Objective gives the total derivative of the objective with respect to a declared parameter,

df/dp  =  df/dp|_x  +  sum_i (df/dx_i)(dx_i/dp)

which is the quantity an outer-loop optimization, a design-of-experiments score, or a “which parameter is my objective most exposed to” question actually wants. It is one number per parameter, on the same convention as the rest of the call: pass the Pyomo object, get a float.

Both halves are included. A parameter that appears in the objective contributes its explicit partial as well as its effect through the solution — on min (x - p)^2 + 3 p^2 subject to x + y == 5, where the optimum sits at x = p, the whole answer is the explicit partial and a chain-rule-only reading would return 0 instead of 6p. Nothing about that 0 looks wrong, which is why pyomo-pounce/tests/test_issue_878_objective_total_derivative.py carries a fixture whose implicit half vanishes.

This works because declare_sens_param has already rewritten the parameter into a variable pinned by a defining equality, so p is an ordinary coordinate of the solve: the objective gradient carries df/dp|_x in p’s own slot and the derivative column carries dp/dp = 1 there. One contraction picks up both terms, with no second index convention to get wrong.

Only the active objective of the solved model is accepted; a deactivated one left on the model from another formulation is refused by name rather than answered with the solved objective’s gradient.

sens_jacobian returns exact first-order derivatives (unit-perturbation backsolves, no finite differencing); sens_solution combines the stored derivative columns for arbitrary perturbed values after the fact. Its perturbation is measured from the solve point, not the Param’s current value, so writing a measurement into the Param before asking (the receding-horizon pattern) does not change the answer. It also warns when the linear step leaves the variable bounds, and mode="fix_relax" pins those variables and re-solves instead, covered in Bending the estimate around a bound below. mode="path" applies the change a little at a time and records where the active set changes along it, covered in Applying the change a little at a time below. There is one exception to the warning, a bound written on a declared Param, covered in Declared Params in variable bounds below. sens_solution_report() measures the same step and reports where the active set changes along it, covered in What the step did about the bounds. Multiplier sensitivities are available for equality constraints. Models without declarations solve through the ordinary AMPL/CLI path, unchanged. See python/notebooks/25_pyomo_sensitivity.ipynb for a worked optimal-control example (initial conditions as parameters; the first-move gradient IS the NMPC feedback gain).

Bending the estimate around a bound: mode="fix_relax"

sens_solution() takes the linear step, and where that step leaves a variable’s bound it clamps the value and warns. Clamping is all the linear step can do, and it costs more than the one variable: every other variable keeps the value the step gave it, computed on the assumption that the clamped one was free to move where the step said. The result satisfies the bounds and no longer satisfies the constraints.

mode="fix_relax" repairs the active set the step implies instead, which is upstream sIPOPT’s strategy of that name and both of its cases. A variable the step carries past a bound is pinned there, activating it. A bound multiplier the step drives negative is set to zero, deactivating that bound so the variable can move. Each adds a row to the held factorization and re-solves, so the other variables move with it:

sens_solution(m, [(m.setpoint, 3.0)])                      # clamps
sens_solution(m, [(m.setpoint, 3.0)], mode="fix_relax")    # pins and re-solves

Both halves matter and they fail differently. On a model where y = 2x + 1 and x hits its lower bound, the linear step returns y = -5, which does not satisfy the constraint at all, while pinning x returns y = 1, matching a full re-solve. On a model whose bound wants to release, the linear step is stuck at x = 0 where the answer is x = 1.667, because the step preserves complementarity and nothing but the release lets the variable off its bound.

Both modes also carry a correction for the barrier. The step is taken against a factorization held at the solve’s final mu, so on its own it estimates where the BARRIER problem’s solution moves rather than the original problem’s, and the two differ by O(mu). That is invisible at a converged tolerance and is not at a loose one: against sIPOPT the uncorrected step differs by 9e-6 at tol = 1e-3 and by 2e-9 at tol = 1e-8. There is no option for it, since there is no reason to want the barrier problem’s answer.

A pass takes every crossing it can see, pins them together and re-solves, which is upstream’s own loop: one violation list per pass, and while (bounds_violated) as the termination condition. Each pass rebuilds the Schur complement over the pins so far, so a pass carrying k of them costs one dense k × k solve and k + 1 back-solves. A pin never rebuilds the factorization, which is what keeps it cheaper than re-solving.

bound_eps sets how far outside a variable bound a step has to end to count as having left it, and so decides what a pass pins, what sens_solution() clamps, and what crossed reports. It is absolute, as the refinement’s own test is. Unset, it is how far outside the solve itself was willing to settle, so nothing moves for a caller who does not set it. A constraint row keeps its own floor, and a bound is released when the step drives its multiplier negative past the solve’s own margin, whatever bound_eps is. mode="path" reads no such margin, and passing it under linear or path warns.

max_pdpert refuses rather than answering when the converged factor carries an inertia correction above the value given, since every sensitivity output inverts that factor and a perturbed one answers for a nearby problem. sens_solution_report().perturbations reports the same numbers for a caller who would rather read them.

predictor_iter caps the passes, and is a safety limit rather than a budget. It was a budget while a pass took only the worst crossing, which needed as many passes as there were crossings — and on a model with more crossings than passes the limit, not the violations, decided where the loop stopped. On the CSTR of notebook 36 that put the pin count at exactly the budget for every budget tried, and at 100 pins (half that problem’s degrees of freedom) the refined step came back 8.6 times worse than the unrefined one (gh#732). sens_solution()’s warning now names which stopping condition was reached rather than inferring it from the pin count.

A release does re-factor, once per released set. It has to. An active bound contributes sigma = z / s to the KKT’s x diagonal, and the tighter the solve the larger that term, so the released system’s information is destroyed in the converged factor to about eps · sigma. Computing a release from the held factorization therefore gets worse the better the solve converged — at tol = 1e-10 the released answer was off by 2e-4 while at a looser 1e-6 it was off by 7e-9. Dropping the bound’s sigma and re-factoring removes the dependence entirely. One factorization still sits an order of magnitude under the twenty to a hundred a re-solve runs, and a step that releases nothing pays nothing. This is the one place the loop departs from upstream, which puts the multiplier’s row in the same violation list as the primal crossings and takes a Schur row over it: that is the computation the eps · sigma cancellation above is measuring. The pins survive a release either way — their right-hand sides are re-measured against the re-solved base rather than the pin set being cleared.

Releasing is the half that has to be careful about how much it does at once, and the asymmetry is worth stating. A pin ADDS a condition, so asking for too many shows up honestly as an augmented system that cannot be solved. A release REMOVES one: each bound taken out of the active set is stiffness that is no longer holding its variable there. Take too many at once and variables that were sitting on their bounds are carried off them, with no degrees of freedom left to pin them back — and no failed solve anywhere to say so. So a release batch is kept only when the step it produces is no further outside the bounds than the step in hand; otherwise the most negative multiplier goes alone and the next pass re-measures the rest under it.

Three things stop it short of holding every bound. The pass limit, which a caller can raise. The problem’s degrees of freedom, which no limit helps: pinning uses one degree of freedom each, and past that no step holds every bound at once, so the pin is refused rather than returned from a singular system. And the refinement ending further outside the bounds than the step it started from, which returns the unrefined step instead — repairing an active set has to beat not repairing it. In each case sens_solution() warns, names the variables still outside, and says which of the three it was. clamp then decides what happens to them, exactly as under linear.

A pass is also refused when its correction is out of scale with the step it corrects, not only when a pinned row misses its target. Checking the pinned rows alone is what let gh#732’s hundred pins each land within 1e-3 of where they were asked to go while the step as a whole came back unusable: hitting the pinned coordinates says nothing about what the correction did to the other thirteen hundred.

What counts as outside a bound is not a tolerance you pass. It comes from the solve, which was willing to leave a converged point bound_relax_factor outside its bound, so anything within that is on the bound rather than past it.

This is what sens_boundcheck turns on for the CLI and the Rust API, and it mirrors upstream sIPOPT’s option of that name.

Applying the change a little at a time: mode="path"

mode="fix_relax" decides every active-set change from full steps taken at the base point. mode="path" follows the solution along the perturbation instead: it takes the largest fraction of the change the current active set allows, applies the one change that happens there, and continues under the updated set. The prediction is piecewise linear in the parameter. For a QP that is the exact solution path, since a QP’s solution is piecewise affine in the parameter. For an NLP the one error left is the linearization at the base point, because nothing is re-linearized between breakpoints.

Three kinds of breakpoint end a segment. A free variable reaches a bound and is held there. A bound active at the base has its multiplier fall to zero and the variable leaves it. A bound the path itself started holding stops binding under a later direction and the variable leaves it again. That last kind is what no decision at the base point can represent: a variable can arrive at a bound partway through the change and depart before the end.

sens_active_set_changes() returns that record, which is the part no other mode produces. It takes the same perturbation argument sens_solution() takes:

from pyomo_pounce import sens_active_set_changes, sens_solution

sens_solution(m, [(m.setpoint, 3.0)], mode="path")
for c in sens_active_set_changes(m, [(m.setpoint, 3.0)]):
    print(c.fraction, c.var.name, c.bound, c.action)

Each entry holds the fraction of the perturbation at which the change happens, the variable, which bound ("lower" or "upper"), and whether the variable "reaches" it or "leaves" it. The first entry’s fraction is how much of the perturbation the held solve’s active set survives unchanged.

Where the two modes settle the same active set they give the same prediction. Where the changes are spread out along the perturbation they differ: on the notebook’s CSTR at a change large enough to release thirteen bounds, fix_relax decides all thirteen at once from base-point multipliers and its prediction lands below even mode="linear" (worst relative miss 0.950 against 0.833), while mode="path" applies each release at the fraction the record names and stays the most accurate of the three (0.626). At changes this large every first-order prediction degrades: the CSTR trajectories read high near the start of the horizon in every mode, which is the base-point linearization and not something more segments repair.

predictor_iter is the same knob it is under fix_relax: it caps the active-set changes applied, and past the cap the rest of the perturbation is taken in one step under the active set reached, with the warning naming the cap. On the cost side a reach adds a Schur row without re-factoring, each release re-factors once, and the wall time grows about linearly with the changes applied, well under a re-solve.

See python/notebooks/36_active_set_parametric_sensitivity.ipynb for the worked CSTR example behind those numbers, including predictor_iter sweeps of both modes against re-solve wall time.

A held solve at a kink: degeneracy

A solve can converge with a bound weakly active, the variable on the bound with a multiplier of the same order as the slack, both of order the square root of the barrier parameter. The solution as a function of the parameter has a kink there, with a different one-sided derivative on each side, and no single linear step is right for both. The activity classifier reports such a bound as weakly_active, or as ambiguous where the coordinate is coupled to a neighbour — see ambiguous is not “probably not a kink”. Everything downstream treats both as weak, so the distinction does not reach the step.

The factorization carries every bound as sigma = z / s on the variable’s diagonal. At a strongly active bound that is around 1e8 and the variable cannot move, at an inactive one around 1e-8 and the bound imposes nothing, and at a kink it is of the order of the reduced curvature along that coordinate — order one on a well-scaled model — so the bound is only partly enforced, which is wrong for both sides. The thresholds that decide activity elsewhere have no answer at a kink, since the two quantities they compare are the same size.

degeneracy on sens_solution(), sens_solution_report(), and sens_active_set_changes() selects what happens then:

sens_solution(m, [(m.p, 2.5)], degeneracy="directional")   # the default
sens_solution(m, [(m.p, 2.5)], degeneracy="one_sided")     # the thresholds' answer
sens_solution(m, [(m.p, 2.5)], degeneracy="release_all")   # released, undecided

"directional" decides each weakly active bound for the perturbation’s own direction by the directional-derivative QP (the sIPOPT paper’s eq. 14). The weakly active rows are released, removing the order-one sigma, in one factorization that serves the whole decision, and the direction of the released system is computed. Rows it moves toward their bounds are the ones the direction engages, and their pin forces solve a small quadratic program, one variable per engaged row with a nonnegativity bound, whose optimality conditions are eq. 14’s complementarity: each engaged bound either holds with a nonnegative force or releases and moves feasibly. An engaged row is decided only when its bound is at a kink, which the QP reads off its own reduced matrix: the barrier weight times the row’s own diagonal is 1 at an exact kink, at any curvature, coupling, or scaling, and falls as the squared ratio of kink width to slack away from one. A row far below one is dropped and its plain movement stands, since a pin force there would hold the coordinate a full slack from where its bound actually is, and the error of not deciding is bounded by that same slack. A row an equality pins is the limiting case, its diagonal exactly zero, dropped by the same test. The active-set QP engine solves the rest, the decided direction is checked against every weak row, and the engaged set grows until no new row violates.

A row engages only when its movement toward the bound exceeds a noise band, the square root of the barrier parameter relative to the direction’s norm, so nothing in the decision reads the perturbation’s size and the answer is linear in the step. The band is there because a weak bound’s slack and multiplier carry an uncertainty equal to their own size, so a movement below it cannot be resolved against the bound, and deciding it exactly would assert precision the solve does not contain.

All three modes consume the decision: linear takes the QP direction itself, fix_relax takes it as the predictor its refinement iterates from, and path starts with the held rows pinned and the left rows in its base-activity table, so a bound that is genuinely active for the first stretch of the perturbation, which happens when the held solve sits inside the ambiguous band rather than exactly at the kink, releases at the fraction where its multiplier reaches zero rather than at the start. The record then carries that departure at its measured fraction.

"release_all" releases every weakly active bound undecided, at one back-solve and no QP: the step is the all-released direction, and a weak bound the perturbation actually holds comes back as a bound crossing for whatever runs next. fix_relax pins it, path walks it and records a return to the bound along the path rather than a decision at the base point, and linear clamps the crossing coordinate, which repairs that coordinate alone and leaves its neighbors carrying the released coupling. The trade is the decision’s cost against downstream repair, and the cost is deterministic and independent of degeneracy_iter, which makes this the option for a kinked base point too large for the engagement’s budget, where "directional" pays the failed attempt and falls back to one-sided anyway. At an exact kink under mode="linear" the holding side’s answer is the released one until the clamp truncates it, where "directional" decides it correctly, so the accuracy-first choice at small kink counts remains the default.

On a coupled model the repair is only as good as the mode’s reach, and the three differ. Measured on the coupled kink of pyomo-pounce/tests/test_degeneracy.py, holding side, exact answer x = 0, y = 1:

fix_relax pins the crossing and re-solves, so it repairs the neighbour too. path re-holds the weak bound at the fraction the walk finds the direction pressing into it, and the coordinates behind it re-optimize under the hold, so it reaches the same answer. linear clamps the crossing coordinate only, and the neighbour keeps the released coupling – that is the documented trade, and on a coupled model linear is the mode it costs something.

path answered the one-sided 2/7 here until gh#852, which was split out of this option’s own review: step_along_path barred every base-active bound from its reach scan, so a perturbation pressing into a weakly active one walked the variable out of its box with no breakpoint to stop it, and only a downstream clamp put it back – moving the crossing coordinate and nothing coupled to it. The repair landed in the walk itself, which both the decided and the undecided callers go through, so "release_all" inherited it.

"one_sided" takes the single-sided value the thresholds produce, bit-identical to the behavior without the argument. On the CSTR held at the record’s first breakpoint, a 2% step toward the steady state puts the thresholds on the wrong side: linear and path miss by 0.0077 with an empty record where directional puts all three modes at 0.0018, and fix_relax reaches 0.0018 either way because its own release test happens to read the right sign there, a favorable read that directional replaces with a guarantee.

Undecided is not the same as unenforced, though. Whichever side the thresholds lean toward, path under "one_sided" keeps every weakly active bound inside the box: a perturbation that presses into one is a breakpoint like any other, the walk takes the bound back there, and the coordinates coupled to it re-optimize behind the hold. Before that (gh#852) the walk saw no breakpoint at all on such a perturbation, and the variable left its box for a caller’s clamp to put back — which moves the crossing coordinate and nothing else, so on min (x - p)^2 + 0.1 (y - 1)^2 with y = 2x + 1 and x >= 0, held at the kink p = 0, a step to p = -1 came back with x = 0 against a y of 2/7 where the answer is 1. What "one_sided" gives up at a kink is the choice of side, not feasibility; on that model path and fix_relax now both reproduce the re-solve, and linear is the one that still cannot, since a clamp is all it has. The CSTR figures above are unchanged by it: there the thresholds’ bound is one the step leaves, not one it presses into.

The cost is gated by the condition, and budgeted by degeneracy_iter (default 16): the released solve, one further back-solve per engaged row, and one more to recover the direction all count against it, so the decision costs a handful of back-solves at a kink that engages a handful of bounds. A decision whose engaged set grows pays that recovering solve once per pass, so a set reached in two passes costs one more than the same set reached in one.

A direction that engages more rows than the budget covers falls back to the one-sided step with a warning. Only a budget of zero fails before any work: which rows engage is not known until the released solve has run, so a budget too small to finish still pays that one factorization before reporting the shortfall. The warning names the engaged count and the number to raise degeneracy_iter to, which is the retry price and is a floor, since a later pass can engage more rows. It is always strictly above what the failed call spent, so each retry buys progress. predictor_iter keeps its meaning as the mode’s own work and plays no part in the decision. Detection also returns nothing on a solve with relaxed bounds, where the classifier cannot read the slacks.

sens_jacobian() cannot take a side, since it is asked for a derivative without a direction, so at a degenerate base point it warns, names the variables and bounds, and returns the one-sided value. The direction-aware answer is sens_solution()’s.

Refining the step: corrector_iter

Every mode returns a step, and that step leaves a residual in the barrier KKT system at the perturbed parameter values. Newton iterations drive that residual down against an operator assembled at the predicted point: the Hessian, the constraint Jacobians, and the barrier diagonal are all evaluated at the stepped iterate with the step’s own multipliers, one factorization is paid there, and each iteration afterwards costs one back-solve. A chord iteration contracts at the rate the distance between its operator and the true Jacobian sets, and the predicted point is where the truth is. Under a limited-memory solve the quasi-Newton matrix is kept, since no exact Hessian exists to evaluate elsewhere. corrector_iter is how many iterations to run, on sens_solution() and sens_solution_report(), and it stops early when an iteration fails to improve the residual, so it is a budget rather than a count. It defaults to zero.

The correction aims at the barrier solution at the mu the solve finished on, not at a re-solve, so the accuracy it can reach is bounded by that offset. It does not converge to the exact answer and does not claim to.

What lets it work past a bound crossing is that the predictor already decided which bounds moved, and the corrector applies that decision once before iterating. A bound the step takes off its minimum comes out of the operator, its multiplier held at zero and its complementarity row gone. A bound the step brings onto its minimum has its diagonal raised to the stiffness the barrier assigns there. Every other row carries the predicted point’s own term, in the same frame as the rest of the operator. Both directions are the same change to one diagonal, so the single factorization at the predicted point serves the whole correction.

That decision is where the modes start from different places. fix_relax and path compute an active set and hand it over. mode="linear" holds the active set fixed as it builds the step, so all the correction has to work with is whatever the clamp left sitting on a bound. On the CSTR at a quarter of the change to its steady state that is one bound against the seven the other two pass over, which is why the linear estimate stays furthest from a re-solve. Below the first crossing all three are the same step. A release no step endpoint shows is applied by no mode: the correction can move such a variable partway off its bound, on the weak diagonal entry the step’s clamped multiplier builds at the predicted point, and the estimate is then not the re-solve.

How far the correction reaches is set by how many crossings the predictor hands over rather than by the size of the perturbation directly. Past the crossings the step decided, what limits it is the multiplier handoff below.

The reason it stops is the multipliers rather than the operator. They arrive extrapolated over the whole perturbation, nothing sets them at handoff, and once the perturbation is large that is the dominant error. Fitting them to minimize the stationarity residual at the predictor’s variables does not help: it absorbs the error into the multipliers and removes the signal the iterations need, which is why the algorithm uses that estimate only to initialize multipliers before its first iteration.

So a budget past the crossing count the correction carries buys little, and at large perturbations it can return an estimate no better than the step it was handed. sens_solution() warns when a correction ends without at least halving the residual, so an uncorrected step is never passed off as a corrected one, and sens_solution_report(corrector_iter=...) carries the iterations spent, the residual before and after, and that residual split into stationarity, feasibility and complementarity. The three carry different units and different consequences: a correction can leave the model’s equations nearly satisfied and the multipliers complementary while the Lagrangian’s gradient is far from zero, and only the first two say whether the values can be acted on.

What the step did about the bounds: sens_solution_report()

The clamp warning names the variables it clamped and stops there. sens_solution_report() takes the same perturbation argument sens_solution() takes and measures the same step, so a caller can see how far along the perturbation the active set changes:

from pyomo_pounce import sens_solution_report

r = sens_solution_report(m, [(m.setpoint, 3.0)])
r.alpha          # fraction of the perturbation that fits before a
                 # bound is reached; inf when none lies in the way
r.first          # which variable or constraint is reached there
r.crossed        # {var data: distance past its bound} for the full
r.crossed_rows   # step, and the same for inequality constraints
r.violation      # constraint violation at the predicted point
r.activity       # per coordinate: inactive / weakly_active /
r.row_activity   # strongly_active / ambiguous / unidentified / ...
r.refined        # {name: (before, after)} for every class the
                 # reduced curvature re-classified

Read refined before acting on a class. The classifier that produces activity normalizes a variable’s barrier diagonal by the Hessian’s diagonal and a row’s by the curvature along the row’s own gradient, while the multiplier that produced it is generated by the reduced curvature. The ratio is reduced/diagonal, which is μ-independent — so a genuine kink whose coordinate is coupled reports "ambiguous" and stays there however tightly the problem is re-solved. Coupled coordinates are routine on a collocation model, and reading that class as “probably not a kink” is the inference that shipped gh#763.

The report re-classifies those entries with the reduced curvature before returning, at one back-solve per ambiguous entry and none when there are none, and refined names each one that moved. Pass refine_activity=False to skip it and take the cheap verdict as-is.

What that costs scales with the ambiguous population, not the model size, and on a collocation model the two are far apart. Measured in review of gh#889 on a 62k-variable Radau collocation column: 675 entries were ambiguous, each costs about 29 ms, and the call runs 0.67 s with the refinement off against 20.2 s with it on. Budget it as ambiguous × one back-solve, and read len(rep.refined) afterwards for what a given model actually spent.

Skipping is not free either, which is the half worth stating: the entries that come back "ambiguous" unrefined are a mixture of genuine kinks and genuine non-kinks that no tolerance separates, so the cheap band is wider in the direction that matters. python/pounce/examples/asnmpc_cstr.py calls the report inside a latency-measured control loop and its online guard is narrower because the refinement runs — there, switching it off widens the guard rather than only making the call quicker. Decide it on what the class is for, not on the timing alone.

refine_stop says why the "fix_relax" refinement stopped, one of "settled", "iteration_limit", "degrees_of_freedom" or "worse_than_plain", and is None under the other two modes. A pass pins every crossing it sees, so the pin count says nothing about which limit was reached and this is the only thing that does.

mode and predictor_iter select which step is measured and match sens_solution()’s arguments of the same names. violation and corrector are properties of the step, so a fix_relax estimate needs sens_solution_report(mode="fix_relax") to be described by its own numbers rather than the linear step’s.

Under "fix_relax" and "path" the step stops at the bound, so alpha is 1.0 and crossed is empty for every model. What those two did about the bounds is what "linear" reports at the same perturbation. activity, row_activity and mu come from the converged base point and do not depend on the mode.

alpha comes from a ratio test along the step. Coordinates already on a bound take no part in it, on that side: the gap left at an active bound is the slack the barrier leaves rather than room to move, so scoring it divides two small quantities and would become the minimum on any model carrying an active bound. Which coordinates those are comes from the same classifier Activity classification describes, and for the ones it declines to rule on, from the size of the gap measured against sqrt(mu).

Three further fields say what separates this prediction from the exact value at the perturbed active set, which is what a caller needs when the estimate and a re-solve disagree: mu, the barrier parameter the factorization sits at; perturbations, the factor’s inertia corrections, non-zero when it was regularized; and bounds_relaxed, true when the solve ran with a non-zero bound_relax_factor. That last case is reported rather than raised on, and it empties the two classification maps, because relaxed bounds shift the slacks the classifier reads.

violation is the primal half of the residual. The dual half needs the multipliers at the perturbed point and belongs to a corrector step, which holds them.

Declared Params in variable bounds

A limit is often most naturally written as a bound rather than a constraint:

m.u_max = pyo.Param(initialize=1.0, mutable=True)
declare_sens_param(m.u_max)
m.u = pyo.Var(m.t, bounds=(0, m.u_max))   # the cap, as a bound

A Param left in a bound would be written to the .nl file as a constant at its pre-perturbation value, so the bound would never move and sens_jacobian(m.u[t], wrt=m.u_max) would read exactly 0.0, a wrong answer indistinguishable from a legitimate insensitivity.

POUNCE rewrites such a bound as a constraint over the substituted variable at declaration, so both spellings of the same limit give the same derivative. Expression bounds work too, e.g. bounds=(0, 2 * m.p + 1). Two kinds of variable are deliberately left alone: fixed Vars, whose bounds the solver never enforces, and Vars on deactivated Blocks.

This is a deliberate divergence from pyomo.contrib.sensitivity_toolbox, which substitutes declared Params in constraint and objective expressions only and so still reports zero for the same model; see Compared with the Pyomo sensitivity toolbox. Four things follow from the divergence:

• The bound is dropped from the Var. m.x.ub reads None after the declaration and the NL row carries the reader’s no-bound sentinel 1e19, which is finite, so an isinf() test will not catch it. The moved bound lives on as a row of the _pounce_sens_defs block, part of the declaration’s in-place rewrite.
• sens_solution() does not clamp against a rewritten bound, and raises no clamp warning for it. That is correct rather than an oversight: the bound now moves with the perturbation, so the linear step already respects it to first order.
• sens_covariance()’s bound-active projection still fires. The value the bound held at the solve point is recorded and read back for the activity test, so a declare_sens_fitted variable capped by a declared Param is still projected and still warns.
• It costs a row. A simple bound is handled directly in the barrier; a general inequality costs a slack and a Jacobian row. A model with many Param-dependent bounds trades roughly one row per bound. Only models that write a bound in terms of a declared Param pay this.

A Param pinned to exactly a bound

A related case, and one that used to be silent. A declared Param can pin a variable through an ordinary equality:

m.zc0 = pyo.Param(initialize=1.0, mutable=True)
declare_sens_param(m.zc0)
m.zc = pyo.Var(m.t, bounds=(0, 1))
m.zc_init = pyo.Constraint(expr=m.zc[0] == m.zc0)

d zc[0]/d zc0 is 1 by construction: the equality is linear and says so. When zc0 sits strictly inside zc’s box that is what comes back. When it sits on a bound — zc0 = 1.0 here, or 0.0 — the variable is held by the bound and the equality at once, the force that holds it has no unique split between them, and the solve lands with a bound multiplier far larger than the geometry needs over a slack near roundoff. The barrier diagonal Σ = z/s is the product of both, and it can reach 1e27 against Jacobian entries of 1.

At that point the constraint rows through the variable stop being representable against its own diagonal, and before #737 the whole derivative column read 0.00000 — sens_solution() returned the baseline value, and nothing warned. Σ is now capped at the stiffness those rows can still be seen against, so the equality is enforced again and the column reads what the model says. The cap is a ceiling and not a release: a bound that genuinely holds a variable still holds it, to within roundoff of the variable’s own scale, and a bound-pinned variable that appears in no constraint row is not capped at all — there is no row for its diagonal to swamp, and there the stiffness is exactly what Crossover and the barrier diagonal wants every digit of.

The ceiling holds on every diagonal the sensitivity system builds, the one corrector_iter assembles when a step brings a variable onto a bound included — a bound the corrector newly pins arrives as mu / s² off the step’s own endpoint, which is the same quantity by another name.

It holds on the way back out, too. Folding a bound row into the diagonal is only half of a solve: the row’s own multiplier is recovered from that diagonal afterwards, and the recovery has to divide by the same stiffness the fold used. Before #828 it divided by the uncapped one, so a capped bound was held softly and read back stiffly, and the returned bound-multiplier derivative came out wrong by the cap’s ratio — 1.8e7 against a true 0 on that issue’s fixture, growing as the row’s Jacobian coefficient shrank. corrector_iter then opened on a stationarity residual of the same size, could not reduce it in a single step, and returned the step it had been handed at every budget: the refinement unavailable in exactly the stiff, tightly bounded regime a caller reaches for it in. The multiplier rows now come back through the same cap, on the returned step and on the corrector’s own operator alike, and where the ceiling does not bind nothing moves.

Nothing about the solve changes; this is the sensitivity system only.

Solver options and warm starts

Solver options reach the in-process path the same two ways they reach an ordinary solve: factory-level (SolverFactory("pounce", options={...}) or solver.options[...]) and per-call (solve(m, options={...})), with the per-call mapping winning on conflict. Everything the CLI accepts works here: tolerances, max_iter, scaling, warm-start knobs.

With warm_start_init_point=yes (Python True works too) among the options, the initial multipliers come from the model’s suffixes, the same ones the ASL path uses: dual for equality multipliers, ipopt_zL_in / ipopt_zU_in for bound multipliers, matched by component name (the declaration’s in-place rewrite keeps every constraint’s name, so suffixes keyed by your own constraints match directly; only a call-time sens_params clone still goes through an internal alias). Sign conventions are handled: dual holds the AMPL marginal and ipopt_zU_in Ipopt’s negative-at-upper value, and both are translated to the solver’s internal conventions on the way in.

One deliberate improvement over the ASL path: entries you do not supply take the solver’s own default initialization rather than zero. Through a dense ASL array an absent entry is indistinguishable from a zero multiplier, and a zero bound multiplier on an active bound is a contradictory KKT certificate the solver must first recover from. A suffix knows which entries exist, so an explicit zero is honored (then floored at warm_start_mult_bound_push, exactly as a round-tripped inactive multiplier is) and absence means “initialize as you normally would”: the solver’s own bound_mult_init_val for bound multipliers, and for equality duals the warm path’s 0, which is not the cold path’s least-squares estimate. Seed everything from a prior solve and the two paths behave identically; seed partially and the in-process path degrades gracefully.

Watching the solve (tee=True)

SolverFactory("pounce").solve(m, tee=True) streams the solver’s full Ipopt-style log — banner, problem statistics, iteration table, and end-of-run summary — live to standard output, including inside a Jupyter notebook cell. The log is emitted by the engine itself (the same blocks the pounce CLI prints), so the in-process path just tails it: a long solve shows its iteration table as it runs rather than as one block at the end. Without tee=True the solve is silent, matching the Pyomo convention.

The convex arm (QpSensitivity)

Everything above drives the NLP filter interior-point solver. An LP, convex QP or conic program has its own sensitivity, QpSensitivity in pounce-convex, and the two share their decision-making core (pounce-sens-core) rather than reimplementing it: the same SensBacksolver trait, the same fix-relax / path / directional machinery, the same activity-classification rule. That sharing is the point — it is what stops the two arms drifting on what a kink is.

What the convex arm does that this page’s machinery also does:

What stays NLP-only, and why each is a capability rather than an oversight:

• The corrector. Its entry points take the concrete PdSensBacksolver and read the filter-IPM’s eight-block compound iterate; nothing on the trait describes that shape.
• The covariance and identifiability statistics below.
• The reduced Hessian. Both arms have one, and they are different computations behind one word: sIPOPT’s Schur route here, a null-space projection there. They are deliberately not unified, and the CLI routes a --compute-red-hessian request to this arm for exactly that reason.

The gh#763 rule holds on both arms, and it is the thing to know before reading any status either produces: AMBIGUOUS is not “probably not a kink.” A genuine kink lands there whenever its coordinate is coupled, because the cheap classifier normalizes by a diagonal (a variable) or by the curvature along the row’s own gradient (a row) while the multiplier is generated by the reduced curvature. The ratio is reduced/diagonal, which is μ-independent — re-solving tighter does not separate it. Solver::reduced_activity / Solver::reduced_row_activity answer it at one back-solve per entry, and pounce.sensitivity.solution_report now spends them automatically on the ambiguous entries, reporting what moved under SolutionReport.refined.

From the CLI, a .nl carrying the sIPOPT suffixes on an LP or convex QP is answered on the convex path rather than rerouted; see LP/QP routing.

Parameter covariance and identifiability

For a parameter-estimation model whose objective is a plain sum of squared residuals, the factorization from ONE ordinary solve yields the asymptotic covariance of the fitted parameters. Declare the fitted variables (they stay free) and the residual container while building the model, solve, and ask:

from pyomo_pounce import (declare_sens_fitted, declare_sens_residual,
                          sens_covariance)

m.A = pyo.Var(); m.k = pyo.Var()        # the fitted parameters, free
declare_sens_fitted(m.A, m.k)

m.r = pyo.Var(m.I)                      # residuals, one per data point
m.res = pyo.Constraint(m.I, rule=...)   # r[i] == y[i] - model(A, k, t[i])
declare_sens_residual(m.r)

m.obj = pyo.Objective(expr=sum(m.r[i]**2 for i in m.I))
pyo.SolverFactory("pounce").solve(m)    # one solve

cov = sens_covariance(m)                     # no further information needed

cov[m.A, m.k]               # covariance entry (either order)
cov.std_err[m.k]            # standard error of one parameter
cov.correlation[m.A, m.k]   # correlation matrix entry
cov.matrix                  # dense numpy array, ordered like cov.params
w, V = cov.eigen()          # eigendecomposition, for identifiability

The recipe: the parameter block of the inverse KKT matrix, one backsolve per parameter against the held factor, equals the inverse reduced Hessian of the eliminated problem, and for a sum-of-squares objective cov = 2 sigma^2 (K^-1)_pp. The factor 2 belongs to the unscaled sum of squares (a Gaussian negative log-likelihood objective, SSR / (2 sigma^2), would drop it). The scaling is pinned by test against the analytical linear-regression covariance sigma^2 inv(X^T X) (pyomo-pounce/tests/test_covariance.py).

The noise variance comes from, in order of precedence: sigma_sq= (known measurement variance); the declared residuals (estimated as SSR / (n - n_params), with both numbers derived from the container); or the n_data= fallback for models without explicit residuals, whose SSR is the objective value at the solve — like sens_solution()’s baseline, writing into the model afterwards (a measurement, a warm start for the next horizon) does not move the answer. The solve warns if the declared residuals do not reproduce the objective value (weights or regularization terms would silently corrupt the estimate).

Groups. declare_sens_residual(m.r_conc, group="conc") partitions residuals into noise groups by arbitrary user strings: containers sharing a group (or all ungrouped containers) pool into one estimated variance; distinct groups get their own (cov.sigma_sq becomes a dict), and the covariance switches to the heteroscedastic sandwich form, whose per-group pieces come from the same backsolves. When groups genuinely differ, weighting the objective itself (dividing each group’s residuals by its sigma) is the statistically efficient fix; the sandwich is the truthful report on the unweighted fit.

cov.eigen() returns ascending eigenvalues and matching eigenvectors. An eigenvalue much larger than the rest flags a poorly identified problem: its eigenvector is the parameter combination the data cannot pin down, and the corresponding cov.correlation entries approach +/-1. The returned signs follow the project-wide eigenvector sign convention — the largest-magnitude component of each eigenvector is positive, ties broken by the earliest position in cov.params — so the direction reproduces across machines instead of coming back as whatever LAPACK’s build chose. sens_information().eigen() is the same. sens_covariance warns when the held factor carries inertia-correction perturbations (typically an exactly unidentifiable parameterization) and when the covariance diagonal comes out negative (not a least-squares minimum).

Bound and constraint activity is classified from the solve’s own barrier geometry, not a slack threshold. A STRONGLY ACTIVE bound pins its parameter: zero variance, correlations 0, conditional on the bound, warned. A WEAKLY ACTIVE bound (slack and multiplier vanish together) is KEPT at its full finite variance, corrected for the barrier weight the held factor carries; AMBIGUOUS (loosely converged) and UNIDENTIFIED (curvature below the model’s own noise scale) stay in the free block, each with a warning. A strongly active inequality CONSTRAINT over the fitted parameters pins a combination rather than a coordinate: the matrix is projected on the constraint’s null space, going singular by one per binding row, and the warning names the constraint, the pinned combination, and its conditional information. The same limit written as a bound or as a row returns the same matrix. A binding row that reaches the fitted parameters through free eliminated variables cannot be represented by a restricted normal and is kept unprojected with an explicit warning.

To classify honestly, the declaration-triggered solve sets bound_relax_factor = 0 (slacks must measure distance to your own bounds). This applies to every solve routed through the sensitivity session, not only ones that end in sens_covariance(). If you need the relaxation, pass bound_relax_factor explicitly in options=: your value wins, and sens_covariance() then refuses with a clear error rather than classifying against shifted slacks.

The AMBIGUOUS class is the one this machinery cannot argue away: the interior iterate simply does not carry enough information to decide whether the constraint binds. Crossover (crossover=yes) attacks that directly — it pivots to the active-set path after convergence and returns a point at which a linearly independent set of constraints holds with equality, collapsing the ambiguity into a STRONGLY or WEAKLY ACTIVE verdict. It is a different remedy to the same problem the bound_relax_factor = 0 rule above addresses, and the two compose — genuinely independently, since #654; see Crossover and the barrier diagonal for the measurement. A crossed-over point sits on the declared bounds, i.e. bound_relax_factor inside the box the barrier measured against, so its Σ = z/s used to read z/δ and hold the bound more loosely than an interior iterate would have — degrading, rather than improving, every quantity read off the held factor unless the relaxation was also switched off. Σ is now re-measured against the declared bounds whenever crossover is accepted, so the two options no longer have to be set together.

classify_activity() still requires bound_relax_factor = 0, for the separate reason above: the central-path checks it makes read the barrier’s own slacks, which the relaxation shifts. A solve routed through the declaration-triggered path already sets it; a Solver or SensSolve session you configure yourself does not, unless you ask.

Relation to pounce.curve_fit. This uses the same scale-and-invert-the-reduced-Hessian recipe as pounce.curve_fit — both read a reduced-Hessian block from the held KKT factor and scale it by 2 sigma^2 with sigma^2 = SSR / (n - p) — but with one substantive difference for nonlinear models: curve_fit factors the Gauss-Newton Hessian (pcov = 2 sigma^2 (J^T J)^-1, the expected-information / scipy / pycse.nlinfit convention, always positive semidefinite), while sens_covariance() here feeds the exact Lagrangian Hessian through the .nl bridge, so it reports the observed-information covariance — the full reduced Hessian including the residual-curvature term that Gauss-Newton drops. The two are identical for linear models and in the small-residual / large-n limit, and differ by O(residual x model curvature) otherwise (a few percent on a strongly-curved fit). Neither is uniquely “correct”: Gauss-Newton is the conventional, robust default (it cannot produce a negative variance); observed information is the honest local curvature of the objective you actually solved (Efron & Hinkley 1978) but can go indefinite — which is what the negative-variance warning above is telling you. sens_covariance() offers both: the default hessian="lagrangian" inverts the exact reduced Hessian of the Lagrangian, and sens_covariance(m, hessian="gauss-newton") rebuilds the expected-information form from the residual Jacobian, recovered from the same backsolves at no extra solve (declared residuals required). Reach for it when the numbers must match scipy/nls, when sens_covariance() warns about a negative diagonal, or when the covariance must stay positive semidefinite by construction, e.g. feeding an arrival-cost update in moving horizon estimation. The other difference is the input surface. curve_fit(f, xdata, ydata, ...) is the batteries-included fitter for a callable model f(x, *params) and data arrays: it chooses a starting point, offers robust losses, per-point sigma weights, confidence intervals, prediction bands, dpopt/ddata, and out-of-core streaming, and it projects the covariance onto the active-constraint nullspace when a parameter sits on a bound. sens_covariance() is the post-solve primitive for a model you have already written in Pyomo — residuals as constraints, arbitrary surrounding structure — where you want the covariance of the fit as posed without re-expressing it as f(x, *params). Use curve_fit when the fit is naturally a model-plus-data call; use sens_covariance() to interrogate an existing Pyomo estimation model. Both project a bound-active fitted parameter onto the active-constraint nullspace: sens_covariance() reports the covariance conditional on the active bound (zero variance in the pinned direction, computed by inverting the free block of the information matrix) and still warns, since boundary asymptotics are nonstandard. Only variable bounds on the fitted parameters themselves are detected; a parameter held at the same value by an active constraint row is treated as free (#362). A bound rewritten into a constraint by the rule in Declared Params in variable bounds is the one exception: the value it held at the solve point is recorded, so it is still detected and still projected.

Relation to pyomo.contrib.parmest. parmest is an estimation workflow: multi-experiment data management, bootstrap resampling, and likelihood-ratio confidence regions, at the price of restructuring the problem into its experiment framework, with covariance computed by finite differences or an ipopt re-solve. sens_covariance() is a post-solve primitive: the model as written, one declaration per component, the asymptotic covariance and identifiability diagnostics from the factorization the solve already produced. Use parmest for multi-experiment campaigns and non-asymptotic intervals; use this to interrogate the fit you already have.

See python/notebooks/26_parameter_covariance.ipynb for a worked example with a Monte Carlo validated confidence ellipse and an identifiability diagnosis.

Activity classification

Which bounds and constraint rows are actually holding the solution is a question the converged iterate answers only ambiguously: at a weakly active bound the slack and its multiplier are both O(√μ), so no fixed threshold on either one alone separates “just touching” from “not binding”. Solver.classify_activity() keys on the ratio of barrier curvature to the model’s own curvature instead, which is O(μ), O(1) and O(1/μ) in the three regimes:

solver = pounce.Solver(problem)   # problem.add_option("bound_relax_factor", 0.0)
x, info = solver.solve(x0=x0)

rep = solver.classify_activity()
rep["var_status"]        # ["inactive", "unbounded", "fixed", "strongly_active"]
rep["row_status"]        # ["equality", "strongly_active"]
rep["var_ratio"]         # the ratio behind each call (NaN where nothing was classified)
rep["mu"]                # the barrier parameter the calls were made at

Statuses are inactive, weakly_active, strongly_active, ambiguous (the ratio fell in a gap where this μ cannot decide — re-solve tighter, but see ambiguous is not “probably not a kink” for the case a tighter solve does not separate), and unidentified (the curvature is below noise scale, so the question does not arise). unbounded, fixed and equality mark entries with no barrier geometry to classify.

Both arrays are indexed in user space: var_* follows your n variables and row_* your m constraints, in your order. A variable that fixed_variable_treatment = make_parameter removed from the solve (lb == ub) reports fixed at its own index rather than shifting everything after it.

Two per-entry flags report on the assumptions rather than the geometry: off_central_path (s·z differs from μ by more than 10× on some side) and contaminated (classified inactive yet carrying barrier curvature well above the O(μ) an inactive bound should have — typically a bound that sits close enough to the optimum to bend it).

Inequality rows classify through the same rule, via the curvature along the constraint normal — which is not a reduced curvature either, so ambiguous means the same not-necessarily-a-kink thing there; see the same holds for a constraint row. That is the point of classifying rows at all: move a bound off a variable and onto a row and the activity disappears from the bound-multiplier view entirely, while any covariance or identifiability heuristic keyed on z alone silently stops seeing it (#362).

The call requires the solve to have run with bound_relax_factor=0 (the Ipopt default is 1e-8) and raises ValueError otherwise: relaxed bounds shift the very slacks the classifier reads. The guard tests the value that solve ran under, so setting the option after the fact does not change the answer — set it on the Problem and solve again.

ambiguous is not “probably not a kink”

For a variable the ratio’s denominator is the Hessian diagonal H_ii. At a kink the multiplier is not generated by the diagonal: it is generated by the curvature reduced along that coordinate, i.e. what is left once the other free variables re-optimize. Eliminating a free partner y from [[h, c], [c, m]] leaves h - c²/m, and Σ equals exactly that, so the ratio is

r = reduced / diagonal

which is 1 only where the coordinate is decoupled. Couple it and a genuine kink falls out of the band and reads ambiguous — at any tolerance, because that r does not move with μ. Re-solving tighter reports the same thing. On a collocation model, coupling between neighbouring coordinates is the normal case rather than a corner (#763).

So do not read the class as an answer to “is this bound sitting at a kink”. reduced_activity() answers that, normalizing by the reduced curvature instead — one back-solve against the held factor per index, so call it over the entries in question rather than over every bounded variable:

rep = solver.classify_activity()
ask = [i for i, st in enumerate(rep["var_status"]) if st == "ambiguous"]
red = solver.reduced_activity(ask)

red["status"]      # ["weakly_active", ...] — the same rule, reduced denominator
red["ratio"]       # Σ/|q_reduced|; 1 at a kink whatever it is coupled to
red["q_reduced"]   # the reduced curvature itself, natural units
red["var"]         # the user variable index each entry answers about

The default stays the diagonal because the reduced normalizer is the reciprocal diagonal of an inverse, and there is no diagonal-of-the-inverse shortcut: classifying every bounded variable that way costs n back-solves, which on a 62k-variable model is no longer a post-solve diagnostic.

The same holds for a constraint row

A row’s ratio does not divide by a diagonal — it divides by the curvature along the row’s own gradient, |∇dᵀH∇d| / ‖∇d‖². That is a genuine directional curvature, strictly better than a bare H_ii, which is why it was not the one #763 fixed. But it is not a reduced curvature either: the other free coordinates still re-optimize, and what is left after they do is what generates the row’s multiplier. So a row’s ratio is reduced / directional, 1 only where the row’s direction is decoupled from the remaining free space, and a coupled row kink reads ambiguous at any tolerance for the same μ-independent reason (#804).

reduced_row_activity() is the row half of the answer, same shape and same cost — one back-solve per row, so call it over the rows in question:

rep = solver.classify_activity()
ask = [j for j, st in enumerate(rep["row_status"]) if st == "ambiguous"]
red = solver.reduced_row_activity(ask)

red["status"]      # ["weakly_active", ...] — the same rule, reduced denominator
red["ratio"]       # Σ‖∇d‖²/|q_reduced|; 1 at a kink whatever it is coupled to
red["q_reduced"]   # the reduced curvature along the UNIT normal, natural units
red["row"]         # the user constraint index each entry answers about

The row’s own value is a coordinate of the KKT system — the slack the barrier acts on, tied to the model by dⱼ(x) = sⱼ — so the back-solve is the same one the variable accessor makes, one block over: a unit right-hand side in the s block, 1/(K⁻¹)_{sⱼsⱼ} - Σⱼ, then ·‖∇dⱼ‖² to put the answer along the unit normal where classify_activity()’s q lives. Equality rows report equality, as in the report; there is no slack and no barrier multiplier pair to classify.

weakly_active_bounds() and everything built on it (the degeneracy warnings, the directional step) already treat ambiguous as weak for exactly this reason, so the mislabeling is not a wrong answer in the step path: it is a misleading answer to a direct question.

The information matrix

sens_information(model) is the un-inverted sibling of sens_covariance(): the reduced Hessian over the declared fitted block, from the same single solve, in natural units with no sigma^2 anywhere. For a homoscedastic Lagrangian fit, sens_covariance() equals 2*sigma^2*inv(sens_information()) on the free block. hessian= selects the observed ("lagrangian", default) or expected ("gauss-newton") form exactly as in sens_covariance().

The Lagrangian form is built by tangent recovery against the held factorization rather than by inverting the covariance back: the K-inverse columns’ x-blocks are T*M, so T = Zx*inv(M) exactly and R = T'HT with the exact Lagrangian Hessian. The barrier weight cancels multiplicatively, so equality and variable-bound activity carries machine precision at any barrier parameter, including on pinned parameters where a subtract-the-barrier route loses log10(Sigma/q) digits. A binding inequality row is the one exception: it couples through its slack barrier and leaves ~1e-6 relative residue at practical barrier parameters.

Membership and warnings follow sens_covariance(). One disposition is opposite by design: a strongly active (pinned) parameter’s entry is S, the reduction onto the pinned set, NOT a zero row — zero information is the opposite of what a pinned parameter carries — conditional on the rest of the pinned set, with zero cross blocks to the free parameters. Binding constraint rows project the free block on both sides (the pseudo-inverse of the projected covariance). An indefinite Lagrangian block is returned as computed with a warning naming Gauss-Newton as the PSD alternative: refusing would withhold the finding that the point is not a minimum or the model is over-parameterized. eigen() reads identifiability directly: a near-zero eigenvalue is a direction the data does not inform; its eigenvector’s sign follows the project-wide convention.

Choosing the block: of=

Both accessors take of= to reduce onto any block of the solve’s variables off the held factor, post-solve; the declared fitted block is the default, so omitting it is exactly the prior behavior. Accepted forms: a Var (scalar or indexed, every member), an indexed slice (m.x[2, :]), a (Var, iterable) pair, data objects, or a list mixing these.

cov = sens_covariance(m)                      # the fitted block, as before
cov_a = sens_covariance(m, of=[m.a])          # one parameter's marginal
band = sens_covariance(m, of=m.r)             # a predicted trajectory
info_a = sens_information(m, of=[m.a])

Each call re-reduces onto its own argument, so one solve serves as many blocks as are asked about, and each block gets its MARGINAL: everything outside it is profiled out, not held fixed. Sigma estimation always divides by the fit’s own degrees of freedom (a property of the solve, not of the question being asked), so a sub-block’s numbers agree exactly with the corresponding entries of the default answer.

A rank-deficient block, one with more coordinates than the fit has degrees of freedom or with linearly dependent coordinates (a duplicated design point), is the trajectory-band case: sens_covariance() returns its (rank-deficient) marginal, 2 sigma^2 M, the confidence band on the fitted trajectory (add the observation noise for a prediction band), with the membership handling bypassed, and sens_information() raises an error pointing to sens_covariance(), since such a block carries no information matrix. For sens_information(), a block that parameterizes the constraint manifold (size equal to the degrees of freedom) gets the exact tangent construction; a sub-block of the fitted set gets its marginal as a Schur complement of the exact tangent R over the fitted block (never inverting a covariance, so a pinned member costs no digits); other blocks reduce off the held factor with the item-1 corrections, which is benign for free coordinates.

One exception is returned rather than hidden: a strongly active variable OUTSIDE the block is not deleted from the factor, so the block’s numbers are the values conditional on that bound, not the marginal over it. The result carries the list as .conditioned_on (empty when there is none); inside-block activity is membership, not conditioning, and is handled as before. The list is decided by the same classification the block members get, applied per candidate as a singleton block, so it is scale-invariant; only near-bound variables pay the extra backsolve.

Keeping and releasing the factor: sens_retain_kkt(), sens_release_kkt()

The solve factors the KKT matrix to solve the NLP; the only question is whether that factor is kept for post-solve queries. Any declaration keeps it. sens_retain_kkt(model) keeps it with no declaration at all, which is what of= queries with nothing declared need: the MHE case, where the arrival state and the parameters are each queried by of= and neither is THE fitted set. It defaults off, so a solve with no sensitivity pays nothing.

sens_retain_kkt(m)
SolverFactory("pounce").solve(m)
arrival = sens_covariance(m, sigma_sq=s2, of=m.x[:, t0])
params = sens_information(m, of=[m.k1, m.k2])
sens_release_kkt(m)          # done asking: give the memory back now

The retention policy in one place: the factor is kept if anything is declared or sens_retain_kkt() was called, and a Covariance or Information result whose lazy conditioned_on has not been read keeps the session alive through its pending computation until first access. sens_release_kkt(model) is the exit: it drops the model’s hold on the factor immediately, freeing the memory, while declarations and the retain flag still apply to the next solve. Release drops the model’s hold, not a result’s: a Covariance or Information with a pending conditioned_on, and a Jacobian (which reads the factor on every lookup), each hold their own reference, so they keep working across the release and keep the factor in memory until they are discarded. Noise is a separate question: sens_retain_kkt() keeps the factor, not a noise model, and with nothing declared fitted the degrees of freedom for a noise ESTIMATE are unknown, so sens_covariance() under retain-only needs sigma_sq=; the estimation routes (declared residuals, n_data=) raise an error saying so.

Like any declaration, sens_retain_kkt() routes the solve through the in-process sensitivity path, whose solve() surface is not keyword-identical to the ordinary subprocess path (for example, load_solutions=False is not honored there). Adding it to an existing script changes how the solve runs, not just what is kept.

Units and NLP scaling

All sensitivity outputs are in natural (unscaled) units. The IPM holds its converged KKT factor in an internally scaled space whenever NLP scaling is active (the default nlp_scaling_method = "gradient-based" fires when an objective gradient or constraint row exceeds nlp_scaling_max_gradient = 100 at the starting point); pounce undoes that scaling in every held-factor back-solve, so dx, kkt_solve, and the reduced Hessian are independent of how the problem was scaled internally (#128).

That covers user scaling too, on all three of its axes. A per-variable scaling_factor is applied as a change of variables x̃ = d ⊙ x below the algorithm, so the held factor is the scaled problem’s; the factors are carried into the same translation, and every accessor answers in your units (#486). The factors a solve ran under are readable back from Solver.nlp_scaling["x_scaling"] (Python) / Solver::variable_scaling (Rust) — diagnostic rather than a correction to apply, since the outputs already carry it.

classify_activity() is scale-invariant for the same reason, and mostly by construction rather than by undoing anything: its ratios are formed so that rescaling a constraint row or the objective leaves them fixed. Writing a constraint as 1000·x ≥ 0 instead of x ≥ 0 does not move a status, and neither does the solver’s own per-row d_scale. A change of variables is the one case the ratios do not absorb on their own — the identification floor is a single number shared across entries, so a non-uniform d would move entries across it — and there the factors are divided out of the geometry before anything is classified, which keeps a status from depending on the conditioning you asked for. The values the report exports follow the natural-units contract like everything else: var_sigma and row_sigma are the barrier diagonals in the model’s own units, row_normal(j) is the constraint gradient with the solver’s per-row scale divided out, and hessian_vec(v) is the exact Lagrangian Hessian times a user-space vector with the objective scale divided out; classification happens on the scaled quantities internally, the report never shows them. reduced_activity() is invariant on the same terms: its q_reduced is a natural-units curvature, and both sides of 1/(K⁻¹)_ii - Σ_i carry the same d²/df, so the subtraction meets in one frame and the ratio comes out where classify_activity()’s does. reduced_row_activity() reaches the same place with more arithmetic in the way: three dg factors meet in one ratio — the exported Σ carries dg², the back-solved (K⁻¹)_{ss} carries dg⁻² through the natural-units conjugation, and ‖∇d‖² is gathered in a frame that still has dg in it — and the row-scaling leg in tests/reduced_row_activity.rs sweeps six decades of dg to pin it.

Variable indices are user-space, factor rows are not. Everything the sensitivity API reports or accepts — the .col file’s order, the activity report’s var_* arrays, row_normal(j)’s entries — indexes the variables you wrote. The converged factor does not: a variable whose bounds are equal is removed from the solve (fixed_variable_treatment = make_parameter, the default), so its column is absent and every later variable sits one row earlier. The two orders coincide exactly when the model has no fixed variable, which makes the difference easy to miss. Translate with Solver.primal_rows(indices) — None marks a removed variable — before indexing a kkt_solve or parametric_step_full result, just as multiplier_rows has always been required for the y_c block.

In particular, for a parameter-estimation NLP with the parameters pinned by equality constraints, -inv(info["reduced_hessian"]) is directly the parameter covariance — no per-problem scale factor, no need to set nlp_scaling_method = "none". (Sign convention: over pin constraint rows, B K⁻¹ Bᵀ equals the multiplier sensitivity ∂λ/∂p = −∂²f*/∂p², hence the minus in the covariance recipe.)

For callers that calibrated against the pre-#128 behavior, the solver-space value and the factors that relate the two are exposed:

• Python: info["reduced_hessian_scaled"], info["obj_scaling_factor"], info["pin_g_scaling"]; Solver.reduced_hessian(pins, scaled=True), Solver.kkt_solve(rhs, scaled=True), and the Solver.nlp_scaling dict ({"obj": df, "c_scale": …, "d_scale": …, "x_scaling": …}).
• Rust: SensResult::{reduced_hessian_scaled, obj_scaling_factor, pin_g_scaling}, Solver::{compute_reduced_hessian_scaled, kkt_solve_scaled, nlp_scaling, pin_g_scaling}, and PdSensBacksolver::solve_scaled_space.

The relation is H_scaled[i,j] = df / (dc_i·dc_j) · H[i,j], where df is the objective scaling factor and dc_i the pin rows’ constraint scaling factors.

One caveat: the IPM’s inertia-correction perturbations (δ_x, δ_s, δ_c, δ_d) are added to the factor in scaled space, so on a problem whose final factorization needed regularization (e.g. linearly dependent pin rows) the unscaling maps a slightly different perturbed system per scaling method. The perturbations are reported — info["kkt_perturbations"] / Solver.kkt_perturbations (Python), SensResult::kkt_perturbations / Solver::kkt_perturbations (Rust) — so a covariance workflow can assert they are all zero before trusting -inv(reduced_hessian); on well-posed estimation problems the final factor is unregularized and the invariance is exact.

Closed-loop advanced-step NMPC

The CSTR case in 36_active_set_parametric_sensitivity.ipynb now closes the loop around the held-factor examples above. Its reusable driver lives in pounce.examples.asnmpc_cstr. The driver uses optional Pyomo integrations that are not installed by the base wheel. On Python 3.10 or newer, install them with pip install pounce-solver pyomo-pounce pyomo-cvp==0.7.2 before importing the example module. The notebook runs 30 one-minute samples for nominal, constraint-switching, and model-mismatch campaigns and compares:

• a fresh nonlinear-programming solve after every measurement;
• the stale predicted solution with no measurement correction;
• a clamped linear sensitivity update;
• fix-relax and active-set path updates; and
• the same path update behind a full-point acceptance guard and fallback solve.

Each sensitivity-policy sample follows the same ordering: solve the next horizon at the predicted state in the background, receive the measurement, update from that solve’s held KKT factor, validate the corrected point, apply the first piecewise-constant control, integrate an independent plant, shift the horizon, and prepare the next background solve. The full re-solve baseline deliberately does its solve after measurement, so its solve time is online latency; background solve time is reported separately for advanced-step policies.

The guard checks the scaled measurement displacement, the corrector’s full-point feasibility, stationarity and complementarity, corrector progress, path budget, predicted temperature, and ambiguous manipulated-variable activity. A rejection performs a fresh solve at the measurement and resets the warm start and factorization. This is an example policy, not a safety certificate: the guard does not replace plant-side interlocks, state estimation, robust constraint tightening, or a deadline-aware real-time scheduler.

The notebook reports IAE, ISE, stage cost, control movement, maximum temperature violation, active-set changes, fallback counts and fractions, solver failures and recoveries, and median/p95 online latency. A failed warm-started controller solve is retried once from a cold model and both the failure and recovery are counted; an unrecovered cold solve aborts the campaign rather than silently applying an unverified control. Accepted corrected trajectories supply the state, derivative, and control warm start for the shifted horizon; rejected guarded points do not. Timing rows carry the POUNCE commit, model revision, tolerance, platform, Python version, and whether a warm-up was excluded. Re-run timing on the target controller hardware; notebook wall-clock values are evidence for the recorded machine, not portable deadlines.

The online sensitivity timing includes both the applied sens_solution() and the separate sens_solution_report() replay needed for the residual diagnostics; the current public API computes their predictor and corrector work separately. The active-set event-ledger replay is excluded from that timer and reported as diagnostic work. Because the illustrative path-budget guard reads that ledger, a production end-to-end guard should add its cost unless it obtains the record from the applied update itself. The stress campaign’s initial concentration bias is 4.5 times the configured local trust scale, so it deliberately forces an outside-the-validity-region fallback. After that reset, later corrected points can be accepted; this is not evidence that the guard detects subtle model mismatch. Measured temperature is capped at the controller model’s upper bound because that model pins its initial state there; real over-temperature handling belongs to a plant interlock, not this example.

The final experiment holds the paper-scale 100-interval model at its first active-set breakpoint and steps in both directions. It shows why a derivative at a kink is directional and why the guard treats ambiguous control activity as a reason to re-solve.

Verification

Three of the four entry points above — the AMPL CLI, the Rust library and Python — are verified against upstream sIPOPT 3.14.19’s parametric_cpp golden output to within roughly 6e-9 per component. The Pyomo layer has no golden cross-check of its own: it drives the same core through the Python entry point, and its tests cover the modelling-layer translation rather than the numbers.

The bound refinement is verified on that same example, which crosses a bound under upstream’s own perturbation: against a full re-solve the refinement lands within 6e-9 where clamping the crossing coordinate is off by 0.12. It is also checked on a model with three degrees of freedom, where three coordinates cross at once and all three pins hold, and for the refusal when the pins would exceed the degrees of freedom.

Both halves of fix-relax and the barrier correction are also checked against sIPOPT 3.14.19 itself, driven through pyomo.contrib.sensitivity_toolbox, on cases built to separate them:

Each case is one the other two do not reach. For how the two interfaces compare feature by feature, rather than number by number, see Compared with the Pyomo sensitivity toolbox. The release case returns x = 0 without it where the answer is 1.667, since the linear step preserves complementarity and holds the variable on its bound. The barrier case differs by 9e-6 without the correction at tol = 1e-3, and by 2e-9 at tol = 1e-8, which is why it is only visible where the solve leaves mu loose.

Compared with the Pyomo sensitivity toolbox

Pyomo ships its own parametric sensitivity interface, pyomo.contrib.sensitivity_toolbox. It computes the same Pirnay–Biegler quantity POUNCE does; the differences are in how the computation is reached and in what is built on top of it. Everything below is measured against Pyomo 6.10.0.

What the toolbox is

Four entry points over three backends:

The first two work by model surgery: clone the model, replace each declared Param with a Var, walk every objective and constraint substituting occurrences, deactivate all original constraints and rebuild them on a new block, add one paramConst equality per parameter, stamp eight Suffix objects, write .nl/.row/.col into a temporary directory, shell out to the solver binary, and parse the answer back out of files on disk. The rebuild is unconditional — the upstream source notes it: “Unfortunate that this deactivates and replaces constraints even if they don’t contain the parameters.”

Capability comparison

Where the toolbox is the right tool

Three cases, and they are real:

• Any Ipopt build. The toolbox is solver-agnostic; pyomo_pounce’s sensitivity requires POUNCE to be the solver. If you are committed to a particular Ipopt/HSL configuration, that decides it.
• No solver at all. pynumero.get_dsdp_dfdp is a pure equality-Jacobian solve for a square system with as many parameters as degrees of freedom. It ignores inequalities, bounds and multipliers entirely — a limitation for an NLP, and the point for a flowsheet you only want to differentiate.
• Incumbency. pyomo.contrib.parmest and pyomo.contrib.doe both call get_dsdp directly.