API Reference

Problem definition

OptimaSolver.OptimaProblem Type

OptimaProblem{T, F, G}

Gibbs-energy minimization problem in the form:

minimize    f(n, p)            (e.g. G(n) = nᵀ μ(n,p))
subject to  A n = b            (mass conservation, m × ns)
            n ≥ ε              (positivity)

Fields

• A: conservation matrix (m × ns), typically integer-valued

• b: RHS vector (m,)

• f: objective function (n, p) -> scalar

• g!: in-place gradient (grad, n, p) -> nothing (∂f/∂n)

• ns: number of species

• m: number of conservation equations

• lb: lower bounds on n (default: fill(ε, ns))

• ub: upper bounds on n (default: fill(Inf, ns))

• p: parameter tuple passed through to f and g!

Element type

T is promoted over A, b, lb, ub and the numeric leaves of p. The last one is what makes forward-mode differentiation with respect to a parameter work: ForwardDiff seeds p with Dual numbers, and if T were inferred from the constraint data alone it would stay Float64, so the solver's gradient buffer would be a Vector{Float64} and the user's g! would fail trying to write a Dual into it.

OptimaSolver.OptimaOptions Type

OptimaOptions

Solver hyperparameters.

Fields

• tol: KKT residual tolerance (default 1e-10)

• max_iter: maximum Newton iterations (default 300)

• warm_start: reuse previous (n, y) as initial guess (default true)

• barrier_init: initial log-barrier weight μ₀ (default 1e-4)

• barrier_min: minimum barrier weight (default 1e-14)

• barrier_decay: barrier reduction factor per outer iteration (default 0.1)

• ls_alpha: Armijo sufficient-decrease parameter (default 1e-4)

• ls_beta: backtracking contraction factor (default 0.5)

• ls_max_iter: maximum backtracking steps (default 40)

• verbose: print iteration log (default false)

• use_fd_hessian: compute Hessian diagonal via finite differences of ∇f instead of the ideal-solution approximation 1/nᵢ (default false). Enable for problems with pure solid or gas species where the true ∂²f/∂nᵢ² = 0, otherwise the approximation 1/nᵢ causes extremely slow convergence.

use_fd_hessian defaults differently here and on OptimaOptimizer

This struct defaults it to false; the OptimaOptimizer(; …) keyword constructor defaults it to true. Building the same optimizer the two ways therefore selects opposite regimes, and neither is right everywhere.

Measured on chemical equilibria, where the amounts span ten orders of magnitude between the solvent and the trace ions:

	pure water	mixed solid/aqueous
false (analytic 1/nᵢ)	[H⁺]/[OH⁻] = 1.000003 ✓	worst ×6181 ✗
true (finite difference)	[H⁺]/[OH⁻] = 3.78 ✗	worst ×19.8

A Hessian affects the rate of Newton's method, not its limit, so a solver that merely converged slowly would still reach the right answer. It does not — the iteration stalls, identically at 300 and at 200 000 iterations — which places the fault in the interior-point loop, most likely in how a pure phase (whose objective is linear in its amount) is driven to its bound. Until that is resolved, choose this flag deliberately rather than relying on either default.

OptimaSolver.OptimaState Type

OptimaState{T}

Mutable solver state — primal variables n, dual variables y (Lagrange multipliers for A n = b), and the barrier parameter μ.

Warm-starting: pass the converged state from a previous solve as u0 to solve; the solver will initialize (n, y) from it.

OptimaSolver.OptimaResult Type

OptimaResult{T}

Immutable solver output.

Fields

• n: equilibrium mole amounts

• y: Lagrange multipliers at convergence

• iterations: total Newton iterations

• converged: convergence flag

• error_opt: final KKT optimality residual

• error_feas: final feasibility residual, each row scaled by its own budget (see row_scales) — this is the quantity tol is compared against

• error_feas_abs: the same residual unscaled, ‖An − b‖∞, for reporting

Canonicalizer

OptimaSolver.Canonicalizer Type

Canonicalizer{T}

Decomposes the conservation matrix A (m × ns) into canonical form:

A Q = [B  N]    with B = basic block (m × m, full rank)

where Q is a column permutation such that the first m columns are linearly independent. The LU factorization of B is cached and reused for every Newton step without re-factorization.

Fields

• A: original conservation matrix (m × ns)

• Q: column permutation (pivot indices into 1:ns)

• Qinv: inverse permutation

• jb: basic variable indices (length m)

• jn: non-basic variable indices (length ns-m)

• B: basic block A[:, jb] (m × m)

• BLU: LU factorization of B (used for Newton and sensitivity)

• R: R = B⁻¹ N (m × (ns-m)), the reduced-cost matrix

• ns: number of species

• m: number of constraints

• rank_A: effective rank of A

Solver

solve is exported and re-exports SciMLBase.solve, so using OptimaSolver is sufficient. solve! (the in-place variant) is not exported; use the qualified name OptimaSolver.solve!(...) or import OptimaSolver: solve!.

CommonSolve.solve Function

solve(prob, opts; u0, y0) -> OptimaResult

Solve the Gibbs minimization problem prob and return an OptimaResult.

Arguments

• prob: OptimaProblem

• opts: OptimaOptions (keyword; defaults to OptimaOptions())

• u0: initial guess for n (keyword; defaults to b/m spread)

• y0: initial guess for y (keyword; defaults to zeros)

Warm-start

Pass a previous OptimaResult as u0 = prev_result and the solver will initialize from prev_result.n and prev_result.y.

solve(prob, can, opts; u0, y0) -> OptimaResult

Variant that accepts a pre-built Canonicalizer (avoids recomputing QR when prob.A is fixed across many solves, e.g. during a temperature scan).

SciMLBase.solve(opt_prob, alg::OptimaOptimizer; kwargs...) -> OptimizationSolution

Convert a SciML OptimizationProblem to an internal OptimaProblem and solve with the OptimaSolver primal-dual method.

The OptimizationProblem is expected to carry:

• f.f: objective (u, p) -> scalar

• f.grad: in-place gradient (g, u, p) -> nothing (or nothing)

• prob.cons: equality constraints (res, u, p) -> nothing (A u = b encoded as residual)

• prob.lcons, prob.ucons: lower/upper constraint bounds (should be equal for equality)

• prob.lb: lower bounds on u

• prob.u0: initial guess

• prob.p: parameter tuple

Gradient fallback

If f.grad is nothing, a ForwardDiff gradient is constructed automatically.

OptimaSolver.solve! Function

solve!(state, prob, can, opts) -> OptimaState

Run the interior-point loop, mutating state in-place.

Arguments

• state: OptimaState — initial iterate on entry, solution on exit

• prob: OptimaProblem

• can: pre-built Canonicalizer for prob.A

• opts: OptimaOptions

KKT solver and certificate

Newton on the KKT conditions, and the proof that a point is optimal. Derived in Newton on the KKT system in multiplier space.

OptimaSolver.SolutionPhase Type

SolutionPhase(members, j_ref; always_present = false)

A set of variables whose hᵢ depends on the composition of the set — a mixing phase — together with the position within members of its reference.

Every member of such a phase satisfies hᵢ = uᵢ − gᵢ, and for all but the reference that condition is inverted directly. The reference is exempt for a structural reason: hᵢ is a logarithm of a mole fraction, hence bounded above, so hᵢ = uᵢ − gᵢ has no solution for an arbitrary y. Inverting it is not merely slow, it can be infeasible, and an inner loop that included it can never report convergence — which invalidates the outer Jacobian, that Jacobian being derived on the assumption the inner conditions hold exactly. Its stationarity is carried by the outer system, where it fixes the phase's total.

always_present = true exempts the phase from the stability test: the aqueous phase of a wet system is there by hypothesis, a solid solution is not.

OptimaSolver.DualNewtonProblem Type

DualNewtonProblem(A, g, h; phases, idx_bounded, params)

A convex program in the form solved by dual_newton_solve:

Arguments

• A: the m × n equality matrix.

• g: the constant part of the gradient, ∇f = g + h(x).

• h: callback h(x, params) -> Vector, the state-dependent part.

Keywords

• phases: the mixing phases, as SolutionPhase values. Their members are strictly positive while the phase exists, hence parameterized by ln x and recovered from their own stationarity.

• idx_bounded: variables outside every phase. Their hᵢ does not depend on the composition — a pure phase, of unit activity — so they are either at a stationarity of their own or at zero, and an active set decides which.

• params: passed through to h.

• conservation_rows: which rows of A x + Aq q = b the degeneracy criterion of degenerate_components may be applied to. Defaults to all of them, which is right when every row conserves an element or the charge. It must NOT be all of them once a row means something else. A reactivity row nᵢ − Σⱼ νᵢⱼ Δξⱼ = nᵢ(0) has a single positive entry on x and, for a product that starts absent, a zero right-hand side — exactly the shape the criterion reads as "no matter of this component exists". It then pins that row's multiplier to DEGENERATE_POTENTIAL and declares the species dead, so a solid product can never form: measured, the stationarity residual sat at 458 and the reaction extents came out 4.3 times short. "This species starts at zero" and "this component is absent from the system" are different statements, and only the second licenses the criterion.

• always_active: bounded variables that are never dropped from the active set. A variable whose amount is fixed by a linear row is not deciding anything by a sign test: it is determined, and its extra multiplier makes its stationarity satisfiable at whatever amount the row demands. Without this it can never enter, because it starts at zero and the drop rule removes anything below si_tol — which is exactly what happens to the products of a solid-to-solid reaction whose extents pin them.

The two kinds of variable

The distinction is not cosmetic. A pure phase satisfies gᵢ + hᵢ = uᵢ when present and ≥ when absent, and it can be exactly zero. A member of a mixing phase cannot: its activity goes to −∞ as its fraction goes to zero, so it is never exactly absent while the phase exists. The active set for a mixing phase is therefore over the PHASE, and the criterion is a tangent-plane test rather than a sign of a saturation index.

OptimaSolver.DualNewtonOptions Type

DualNewtonOptions(; tol, maxit, max_active_updates, si_tol, verbose)

• tol: tolerance on the KKT residual.

• maxit: Newton iterations per active set.

• max_active_updates: how many times the active set may change.

• si_tol: saturation index above which a variable at its bound is admitted.

OptimaSolver.dual_newton_solve Function

dual_newton_solve(prob, b, x0; opts) -> (; x, y, q, active_phases, active, converged)

Solve prob for the right-hand side b, starting from x0.

x0 supplies a neighborhood, not a feasible point: this is a Newton method, and the intended use is to hand it the answer of an interior-point solve. Verify the result with kkt_certificate; for a convex problem that certificate is a proof of global optimality.

Two active sets

Over the bound-constrained variables, on the sign of uᵢ − (gᵢ + hᵢ): a pure phase is present exactly when that index vanishes, absent when it is negative.

Over the mixing phases, on Michelsen's tangent-plane measure Σᵢ exp(uᵢ − gᵢ) − 1: a solution phase forms when a trial composition of it lies below the tangent plane of the current state. That test is what a mixing phase requires, since its members are never exactly absent while it exists and it has no single saturation index.

Both admit one candidate per round, the most violated, and the sets visited are recorded, which bounds the loop by the number of subsets and therefore terminates. Admitting a batch feeds a cycle in which a variable is admitted, driven negative, dropped and readmitted.

OptimaSolver.kkt_certificate Function

kkt_certificate(prob, x, b; floor = 1e-25) -> (; stationarity, feasibility,
                                                worst_violation, n_interior,
                                                n_forced_zero, optimal)

Check the KKT conditions at x, independently of how it was obtained. For a convex problem they are sufficient, so optimal = true is a proof.

What is checked, and on which variables

A variable is INTERIOR when xᵢ > floor. There the condition is the equality ∇fᵢ + (Aᵀy)ᵢ = 0, and y is obtained from those variables by least squares. Below floor a variable is at its bound, where the condition is the INEQUALITY ∇fᵢ + (Aᵀy)ᵢ ≥ 0.

Getting that split wrong is not a detail: imposing the equality on a variable held at 1e-16 whose stationarity value is e⁻³⁰⁰ misstates hᵢ by 263 units, and the check then reports a residual of 74 for a point solved to 5e-12.

Variables carrying a component whose right-hand side has vanished are excluded from both tests: they are zero by the CONSTRAINT, and the multiplier of a component nobody supplies is determined by nothing.

OptimaSolver.degenerate_components Function

degenerate_components(A, b) -> Vector{Int}

Rows k of A x = b whose right-hand side forces every variable carrying component k to zero.

The criterion is not simply b_k ≈ 0. With x ≥ 0, the row Σᵢ A_{ki} xᵢ = 0 forces xᵢ = 0 for all i with A_{ki} ≠ 0 only when the non-zero entries of the row share a sign: a sum of non-negative terms vanishes only if each vanishes. A row with entries of both signs permits cancellation and forces nothing.

That distinction is not academic. In a chemical system the H+ row carries −1 for OH- and +1 for H+, so b = 0 there is the ordinary state of pure water — declaring it degenerate kills the whole acid–base system and returns pH 7.000 with the solid undissolved.

OptimaSolver.stationarity_capacity Function

stationarity_capacity(prob) -> Int

The number of conservation rows m, which is the largest number of simultaneous stationarity conditions the element potentials can carry.

This is Gibbs' phase rule, in the form the dual system takes. A bound-constrained variable held ACTIVE contributes uᵢ = gᵢ, i.e. aᵢᵀ y = −gᵢ, one linear equation in y; a mole-fraction mixing phase contributes logsumexp(uᵢ − gᵢ) = 0, one more, nonlinear but still a condition on y alone. A phase with a solvent does not, because its reference equation involves the composition. With y ∈ ℝᵐ there is no y satisfying more than m of them, so an active set carrying more cannot support a solution — the Newton residual cannot reach zero for any iterate, and the least-squares step merely spreads the violation over the rows.

That is not a slow case, it is an unsolvable one, and it has to be excluded by construction rather than discovered. Measured on an LC³ equilibrium the active set grew to 15 pure phases and 5 solid solutions — 19 conditions on 12 components — and the solve came back with a stationarity residual of 18 and an element balance of 800 having never had a solution to find.

OptimaSolver.DEGENERATE_POTENTIAL Constant

DEGENERATE_POTENTIAL

Multiplier assigned to a component whose right-hand side has vanished. Any value large enough to drive exp(uᵢ − gᵢ) below the floor for every variable carrying that component will do; the solution does not depend on it.

Sensitivity

OptimaSolver.SensitivityResult Type

SensitivityResult{T}

Sensitivity of the equilibrium composition n* with respect to:

• ∂n_∂b: ∂n*/∂b (ns × m) — response to changes in mass-balance RHS

• ∂n_∂μ0: ∂n*/∂(μ⁰/RT) (ns × ns) — response to standard potentials

These Jacobians can be used directly as the Jacobian of the RHS in an ODE coupling kinetics (rates→ b) and thermodynamics (T-dependent μ⁰ → potentials).

OptimaSolver.sensitivity Function

sensitivity(prob, n, y, h; μ) -> SensitivityResult

Compute the sensitivity matrices ∂n_/∂b and ∂n_/∂(μ⁰/RT) at the converged point (n, y) with Hessian diagonal h.

The KKT Jacobian is:

J = [ H   Aᵀ ]
    [ A    0  ]

We solve J [dn/dc; dy/dc] = -dF/dc for each perturbation direction.

Using the Schur complement (same decomposition as the Newton step):

S dy = rhs_dual         S = A H⁻¹ Aᵀ
dn   = -H⁻¹ (ex_rhs + Aᵀ dy)

Arguments

• prob: OptimaProblem

• n: converged primal iterate (ns,)

• y: converged dual iterate (m,)

• h: Hessian diagonal at (n, y) (ns,)

• μ: barrier parameter at convergence

SciML interface

OptimaSolver.OptimaOptimizer Type

OptimaOptimizer

Drop-in SciML optimizer implementing the OptimaSolver primal-dual interior-point algorithm for Gibbs-energy minimization.

Constructors

OptimaOptimizer(; tol=1e-10, max_iter=300, warm_start=true, verbose=false)
OptimaOptimizer(opts::OptimaOptions)

Fields

• options: OptimaOptions with all algorithm hyperparameters

• _cache: Ref{Union{Nothing, OptimaResult}} — previous solution for warm-start

OptimaSolver.reset_cache! Function

reset_cache!(alg::OptimaOptimizer)

Clear the warm-start cache. Call this when the chemical system changes (new set of species, different A matrix).

Internal components

The following symbols are exported for testing and extension purposes. They are not needed for typical usage.

KKT residual and Hessian

OptimaSolver.KKTResidual Type

KKTResidual{T}

Holds the KKT residual vectors and associated norms for one evaluation.

OptimaSolver.kkt_residual Function

kkt_residual(prob, n, y, grad_f, μ) -> KKTResidual

Compute the KKT residual at (n, y) with barrier weight μ.

• grad_f: gradient ∇f(n) evaluated outside (allows caching)

• μ: log-barrier weight (scalar, T-compatible for AD)

OptimaSolver.hessian_diagonal Function

hessian_diagonal(prob, n, μ, hess_f_diag) -> Vector

Return the diagonal of the barrier-augmented Hessian:

H_diag[i] = hess_f_diag[i] + μ / (n[i] - lb[i])²

where hess_f_diag is the diagonal of ∇²f(n) (caller-provided). For a convex Gibbs function with positive curvature, H_diag > 0 always.

Why this is μ/sᵢ² and not the primal-dual zᵢ/sᵢ

Optima and Ipopt build this term from an iterated bound multiplier, Σᵢ = zᵢ/sᵢ, and that is the one substantive difference between a primal barrier method and a primal-dual one. It was implemented here in full — the z iterate, its own Newton step δz = μ/s − z − Σ δn, its own fraction-to-boundary step length, and the κ_Σ safeguard of Wächter & Biegler (2006) Eq. (16) — and then removed, because it is measurably worse on the problems this package exists for.

The reason is structural, not incidental. A pure phase's Gibbs energy is linear in its amount, so ∂²f/∂nᵢ² = 0 exactly and Σᵢ is not a correction to the curvature — it IS the whole curvature in that direction. μ/sᵢ² is then the exact Hessian of the barrier subproblem actually being solved, and Newton's method on it is exact; zᵢ/sᵢ replaces it with a quantity that only tracks the central path approximately, turning an exact Newton step for the pure phases into an inexact one. A primal-dual method earns its keep where ∇²f supplies the curvature and the multipliers carry information the barrier does not — the opposite regime.

Measured, with everything else identical: on the LC³ clay sweep, μ/sᵢ² certified all five replacement levels and zᵢ/sᵢ failed at 30 % (stationarity 2.07e-10 against 2.91e-11, same composition to four decimals); on the coupled calcite trajectory the two were indistinguishable (one solve short of tolerance, element balance 4.81e-11 mol in both); on the three-species ideal solution they agreed to the last digit, cold and warm-started. No case gained.

OptimaSolver.gibbs_hessian_diag Function

gibbs_hessian_diag(n, p) -> Vector

Diagonal of ∇²G for the ideal/dilute Gibbs function G(n) = nᵀ μ(n,p).

For an ideal solution where μᵢ(n) = μᵢ⁰(T,P)/RT + ln(aᵢ(n)):

• Aqueous solvent (mole fraction): ∂²G/∂nᵢ² = 1/nᵢ - 1/n_aq + 1/n_aq (approximately 1/nᵢ)

• Aqueous solutes (molality): ∂²G/∂nᵢ² ≈ 1/nᵢ

• Pure solids/gases: ∂²G/∂nᵢ² = 0 (or small positive for regularization)

For the general case we use finite-difference or AD. Here we provide the ideal approximation H_diag[i] = 1/nᵢ as a sensible default that is always positive definite.

When using AD via ForwardDiff: do not call this function; instead pass hess_f_diag computed analytically or via forward-mode on the gradient.

Convergence and the barrier schedule

The optimality error and the barrier update are derived in The optimality error, and why the obvious one cannot work.

OptimaSolver.is_converged Function

is_converged(kkt, opts) -> Bool

Return true if the KKT residual is within tolerance.

OptimaSolver.should_reduce_barrier Function

should_reduce_barrier(kkt, μ, opts) -> Bool

Return true if the inner loop has converged sufficiently to reduce μ.

We use a relaxed inner tolerance: 10× tol (or the current μ if larger), so that we tighten the barrier aggressively when far from the solution and gently when close.

OptimaSolver.reduce_barrier Function

reduce_barrier(μ, opts) -> μ_new

Apply the barrier reduction schedule.

μ_new = max(barrier_min, barrier_decay * μ)

Newton step

OptimaSolver.NewtonStep Type

NewtonStep{T}

Workspace for the Schur-complement Newton solver. All matrices and vectors are pre-allocated once and reused across Newton iterations.

Fields:

• S: Schur complement A * diag(1/h) * A' (m × m)

• rhs: RHS ew - A * diag(1/h) * ex for the Schur system (m,)

• dn: primal Newton step δn (ns,)

• dy: dual Newton step δy (m,)

• d: equilibration diagonal sqrt.(diag(S)) (m,)

• AoverH: buffer for A ./ h' (m × ns), shared between the Schur build (S = AoverH * A') and the RHS computation (rhs = ew - AoverH * ex)

OptimaSolver.compute_step! Function

compute_step!(ws, can, h, ex, ew) -> (dn, dy)

Compute the Newton step (dn, dy) by Schur complement elimination.

Arguments

• ws: NewtonStep workspace (mutated in-place)

• can: Canonicalizer (provides A)

• h: Hessian diagonal (ns,), all positive

• ex: optimality residual (ns,)

• ew: feasibility residual (m,)

Returns

(dn, dy) — views into the workspace vectors.

OptimaSolver.clamp_step Function

clamp_step(n, lb, dn; τ=0.995)

Scale the primal step so that n + α dn stays strictly above lb, using the fraction-to-boundary rule with safety factor τ ∈ (0,1).

Returns α ∈ (0, 1].

Line search

OptimaSolver.LineSearchFilter Type

LineSearchFilter{T}

Mutable filter for the line search.

OptimaSolver.line_search Function

line_search(prob, n, y, dn, dy, f_val, grad_f, μ, opts; filter) -> (α, n_new, y_new, f_new)

Backtracking line search with filter acceptance.

Starting from α = α_max (from fraction-to-boundary), tries α, β_α, β²_α, … until the new point (n + α dn, y + α dy) is accepted by the filter or the Armijo condition on feasibility is satisfied.

Returns the accepted step size α and the new iterates.

Variable stability

OptimaSolver.classify_variables Function

classify_variables(n, lb, ex; tol_stable=1e-8) -> (js, ju)

Return indices of stable (js) and unstable (ju) variables.

A variable is stable if (nᵢ - lbᵢ) > tol_stable * max_slack.

OptimaSolver.reduced_step_for_unstable! Function

reduced_step_for_unstable!(dn, ju, n, lb; τ=0.5)

For unstable variables (near lower bound), cap the step to move at most τ * slack toward the bound to prevent crossing.

OptimaSolver.stability_measure Function

stability_measure(n, lb, ex) -> Vector

Compute the stability measure sᵢ = (nᵢ - lbᵢ) * |exᵢ|.

Small sᵢ means variable i is both close to its bound AND has large stationarity residual — likely to be at an active bound at the solution.