Module square_system 

Square-system recovery: when the interior-point solve of a square, constant-objective model fails, solve c(x) = 0 directly.

Why

A model with as many free variables as equality rows, no inequality rows and a constant objective is not an optimization problem. It is a system of nonlinear equations wearing an NLP’s clothes, and every feasible point is a global optimum. That is the class an equation-oriented flowsheet initialization hands a solver (IDAES-style: fix the degrees of freedom, solve the square block), and it is the class where an interior-point failure has the least to do with the model: the barrier, the filter and the Hessian approximation all exist to trade objective against feasibility, and here there is no objective to trade.

Measured on square_flowsheet_resto (32×32, f ≡ 0, feasible by construction). Its Jacobian is nearly rank-deficient even at the solution (σ_min ≈ 3.6e-9 against σ_max ≈ 3e6; the near-null direction is a uniform shift of the eight temperatures, pinned only by exp(-6000/T)), so the solution set is nearly non-isolated. Under hessian_approximation=limited-memory POUNCE on FERAL PR#222 runs to the 3000-iteration cap and every rung of the second-opinion ladder fails too; upstream Ipopt 3.14.19 declares local infeasibility at iteration 50. A Levenberg–Marquardt least-squares solve that uses only J solves it from the same start in 8 iterations. The two figures are from two FERALs and both are named: on FERAL PR#222 at default options, max|c| = 1.9e-10 (recorded when the stage was written); on the published FERAL 0.20 this crate builds against, where rung 4 (feral_increase_quality=no) rescues the leg first by default and the stage is reached only with feral_increase_quality_retry=no, max|c| = 2.874e-10 (scaled 1.710e-10), at max_iter 500 and 3000 alike.

The gate

probe opens this stage only when all of the following hold of the TNLP the ladder was handed:

• every row is an equality (g_l == g_u, finite) and there is at least one;
• the number of free variables — x_l < x_u, i.e. after removing the variables fixed_variable_treatment would remove — equals the row count;
• the objective is constant: the model does not prove ∇f varies (TNLP::derivative_proofs), and ∇f evaluates to exactly zero, every component, at the starting point.

The last test is evaluated, not proved, and that is a deliberate choice. A structural proof would be the stronger gate, but it is not available on the path that matters: presolve wrappers decline derivative_proofs (they change what the rows are), and the C, Python and CasADi surfaces hand POUNCE callbacks, not algebra. What makes an evaluated gate safe is that nothing reported rests on it. The gate only decides whether to spend the attempt. The verdict rests on the certificate checked at the returned point x̂ itself (below), which includes ∇f(x̂) = 0 evaluated there: at a feasible x̂ with ∇f(x̂) = 0, (λ, z) = (0, 0) satisfies stationarity and complementarity exactly, so x̂ is a first-order KKT point by the same definition the IPM’s convergence test uses — whether or not f is constant anywhere else.

The algorithm

A bounded, scaled trust-region Levenberg–Marquardt method on min ½‖c(x)‖² over the free variables:

• Step. Moré’s trust-region form (Moré 1978, MINPACK lmder/lmpar): with column scaling D, the step p̃ = D p solves (ÃᵀÃ + λI) p̃ = −Ãᵀc, à = J D⁻¹, with λ ≥ 0 chosen so ‖p̃‖ ≈ Δ (within 10%) by safeguarded Newton iteration on the secular equation 1/‖p̃(λ)‖ = 1/Δ. The Gauss–Newton step is taken when it already fits.
• Linear algebra. Each λ trial factors the quasi-definite augmented system [[λI, Ãᵀ], [Ã, −I]] with the workspace’s sparse symmetric backend (FERAL) rather than forming JᵀJ, which squares the condition number and densifies. One extra back-solve against the same factor gives (ÃᵀÃ + λI)⁻¹ p̃ for the Newton update. At most MAX_LAMBDA_TRIALS factorizations per iteration.
• Scaling. D_j = max over iterates of ‖J_{:,j}‖₂ (Moré’s non-decreasing column norms), which makes the method invariant to rescaling a variable — the fixture’s variables span 1e-6 to 9e6. The rows are not reweighted: equilibrating them by ‖J_i(x₀)‖∞ was measured on the fixture and lost it from every start tried, while the unweighted merit (MINPACK’s, and scipy’s method='lm') solves it.
• Bounds. The trial point is the projection of x + p onto the box, and the gain ratio is computed for the projected step against the linear model, so a projection that destroys the predicted decrease shrinks Δ — which drives λ up and the step toward the scaled (projected) gradient. Every iterate is inside the box by construction. This is the projected Levenberg–Marquardt idea of Kanzow, Yamashita & Fukushima (2004), in a trust-region rather than a line-search globalization.
• Update. Gain ratio ρ = ared/pred; accept when ρ ≥ 1e-4; halve Δ (to at most the step taken) when ρ < 0.25; set Δ = 2‖D s‖ when ρ > 0.75 or the step was Gauss–Newton. A non-finite residual at the trial point is a rejection with Δ quartered.

Termination and honesty

The LM stops as soon as the iterate passes the acceptance test, or fails on an iteration cap (min(max_iter, 200)), on a vanishing projected gradient of ‖c‖² at a non-root (STATIONARITY_TOL), on Δ collapsing to 1e-15·max(‖Dx‖, 1) (a projected stationary point of ‖c‖² that is not a root — a local least-squares minimum, so the original verdict stands), or on an evaluation / factorization failure.

The acceptance test is the IPM’s own, at zero multipliers: Ipopt’s convergence test asks for the scaled overall error below tol and the unscaled constraint violation below constr_viol_tol, and at (λ, z) = (0, 0) with ∇f = 0 the dual and complementarity terms vanish, leaving max_i s_i|c_i| ≤ tol and max_i |c_i| ≤ constr_viol_tol. s_i is the row scaling the NLP arm computes for the same model, by the same pounce_nlp::orig_ipopt_nlp helpers at the same point (ipm_scaling has the details): none is 1; user-scaling and curvature-based are the TNLP’s g_scaling floored at nlp_scaling_min_value; gradient-based — and equilibration-based, which the NLP arm runs as gradient-based — is the row-max rule over the full Jacobian sampled at the TNLP’s starting point with only fixed variables lifted, applied only when some row exceeds nlp_scaling_max_gradient (or a constraint target gradient is set). The one thing not reproduced is a wrapper the application inserts inside optimize_tnlp that moves the start; that can change s_i, but the unscaled half of the test is independent of s_i, so no choice of s_i admits a point with max|c_i| > constr_viol_tol. On top of that the box is satisfied exactly (projection; no bound_relax_factor widening is borrowed), and ∇f(x̂) must be exactly zero at the returned point. If any check fails the stage reports failure and the caller keeps the original verdict untouched.

Budget and callback

The LM honours max_wall_time and max_cpu_time, as a budget measured from the start of the stage — the per-solve budget each option rung’s optimize_tnlp also gets — and reads the clock once per iterate, after the convergence test. It fires the TNLP’s intermediate_callback at every iterate, the start included, with mode = RegularMode, inf_pr the scaled violation, inf_du = max|∇f|, mu = 0, d_norm the accepted step’s ‖s‖∞ and regularization_size the LM parameter; a false stops it. A stop or an exhausted budget ends the stage without a promotion. No live iterate is installed for the callback (GetIpoptCurrentIterate finds none), because there is no IPM iterate.

What is reported

Solve_Succeeded with x̂, λ = 0 and z_l = z_u = 0, and statistics that describe the LM: iteration_count is its iteration count and, when the application records iteration history, iterations holds one record per iterate (the start included, the IPM’s convention; alpha_primal_char L); final_scaled_objective is f times the objective multiplier the NLP arm would have applied. The post-solve computations that read the IPM’s converged KKT factor (sIPOPT, the reduced Hessian) have nothing to read and are not run; the CLI says so. Those are true KKT multipliers, not placeholders: with ∇f(x̂) = 0 they satisfy ∇f + Jᵀλ − z_l + z_u = 0 exactly. They are unique when J(x̂) is nonsingular and no bound is active; when J is singular (as on the fixture, numerically) the multiplier set is a cone containing zero and zero is its minimum-norm member. They are not the multipliers an IPM would have reported for the same point, because the IPM never reached it.

The root may differ from the one an IPM would have reached: a square system can have several, and the fixture’s nearly non-isolated solution set puts the LM root about 11% (relative, in the temperature block) from the generator’s x*. With a constant objective each is a global optimum.

Structs

FeasibilityTest

The acceptance test, in the NLP arm’s terms — see the module docs.

LmHooks

What the LM consults besides the model: a time budget and a per-iterate observer. LmHooks::default is neither.

LmIterate

One iterate of the LM, as handed to LmHooks::on_iterate.

LmProblem

Everything the LM core needs that is not a callback.

LmResult

What an LM run returns.

RunOptions

How run is to treat its caller’s per-iteration surfaces.

SquareSystemOutcome

What run reports.

SquareSystemPlan

A square, equality-only, constant-objective model, read off the TNLP and ready to hand to run.

Enums

LmStatus

How an LM run ended.

SquareSystemDecline

Why probe did not open the stage.

Constants

LM_ITERATION_CAP

The iteration cap is min(max_iter, this). Measured on the fixture from eight starts (x* scaled by 0.1 … 10 per block): 5 to 143 iterations.

MAX_LAMBDA_TRIALS

Factorizations allowed per LM iteration while solving the secular equation for λ (MINPACK lmpar allows 10 Newton steps; the extra one is the Gauss–Newton trial).

SQUARE_SYSTEM_RUNG_LABEL

The label this stage carries in the second-opinion record (tried, promoted_by, the sweep’s 2nd= column).

STATIONARITY_TOL

MINPACK’s gtol test: stop when every column of Ã is orthogonal to c to within this cosine (the columns are normalized by D, so each |Ã_jᵀc|/‖c‖ lies in [0, 1]), counting only directions the box lets the iterate move in. Set at the noise level on purpose: it exists for the exact case — degenerate_start_infeasible starts where J = 0, and without it the method spends 57 iterations halving Δ at a point no step can leave — and must never fire on the near-null directions of a nearly singular but consistent system like square_flowsheet_resto.

Traits

ResidualSystem

A residual system c: R^n → R^m with a fixed sparse Jacobian pattern.

Functions

bounded_levenberg_marquardt

Solve c(x) = 0 over the box by bounded trust-region LM. See the module docs for the method, the scaling and the termination rules.

default_backend

The backend the stage factors with: FERAL, configured from the same feral_* options as the IPM, and with no summary sink, so the linear-solver summary a caller reads keeps describing the interior-point solves rather than absorbing this stage’s factorizations.

is_enabled

Is the stage switched on? Read as a literal tag so the option-wiring scans see it.

probe

Decide whether the model the ladder was handed is in this stage’s class. Costs one get_nlp_info, one get_bounds_info and — only on a model that is square and equality-only — one starting point, one Jacobian structure query and one ∇f evaluation.

run

Run the LM on a probed model and, if the result passes every check in the module docs, finalize it through tnlp as Solve_Succeeded.