SciML Interface

OptimaOptimizer implements SciMLBase.AbstractOptimizationAlgorithm and is designed as a drop-in replacement for IpoptOptimizer inside ChemistryLab.jl.

Usage with ChemistryLab.jl

The block below is illustrative and is not executed when this documentation is built: ChemistryLab depends on OptimaSolver, not the other way round, so adding it to this package's documentation environment would invert the dependency. The blocks further down that build a SciMLBase.OptimizationProblem directly are executed, and are the ones to trust as an API reference.

using ChemistryLab, OptimaSolver

# Build a ChemicalState as usual, then equilibrate with OptimaOptimizer
state_eq = equilibrate(state0; solver=OptimaOptimizer(tol=1e-10, verbose=false))

OptimaOptimizer handles:

• Variable scaling — each species is scaled by its starting value so that all normalized variables are , critical for convergence when concentrations span multiple decades (e.g. a titration from pH 1 to pH 13).

• Cold-start lifting — absent species (at their lower bound) are raised to a rough element-balance estimate before the first Newton step, avoiding the log-barrier singularity.

• Transparent warm-start caching — each converged result is stored; the next call starts from it automatically.

Constructor options

using OptimaSolver
import SciMLBase

alg = OptimaOptimizer(;
    tol            = 1e-10,  # KKT convergence tolerance
    max_iter       = 300,    # maximum Newton iterations
    warm_start     = true,   # cache and reuse last converged result
    verbose        = false,  # print iteration log
    barrier_init   = 1e-4,   # initial barrier weight μ₀
    barrier_decay  = 0.1,    # μ ← 0.1 μ per outer step
    use_fd_hessian = true,   # correct for mixed solid/aqueous (default true here)
)

OptimaOptimizer(OptimaOptions(1.0e-10, 300, true, 0.0001, 1.0e-14, 0.1, 1.0, 8, 0.0001, 0.5, 40, false, true, true), Base.RefValue{Union{Nothing, OptimaResult}}(nothing))

use_fd_hessian differs from OptimaOptions default

In OptimaOptions the default is false (ideal-solution Hessian  ). In OptimaOptimizer the default is true because real chemical systems typically include pure solid or gas species with zero curvature, for which the ideal approximation would cause extremely slow convergence.

Alternatively, pass a pre-built OptimaOptions:

opts = OptimaOptions(tol=1e-12, verbose=true, use_fd_hessian=true)
alg  = OptimaOptimizer(opts)

OptimaOptimizer(OptimaOptions(1.0e-12, 300, true, 0.0001, 1.0e-14, 0.2, 1.0, 8, 0.0001, 0.5, 40, true, true, true), Base.RefValue{Union{Nothing, OptimaResult}}(nothing))

Warm-start cache management

alg = OptimaOptimizer(tol=1e-10, warm_start=true)

# First call: cold start (no cache yet)
sol1 = SciMLBase.solve(opt_prob_pH5, alg)

# Second call: warm-starts from sol1
sol2 = SciMLBase.solve(opt_prob_pH6, alg)

# Reset cache when the chemical system changes
reset_cache!(alg)

# Next call is a cold start again
sol3 = SciMLBase.solve(opt_prob_new_system, alg)

Non-converged solutions are never cached: if a call fails, the cache is unchanged and the next call falls back to a cold start from opt_prob.u0.

Direct use with SciMLBase.OptimizationProblem

OptimaOptimizer works with any SciMLBase.OptimizationProblem encoding mass-balance constraints as equalities:

ns = 3
A = ones(1, ns)
b = [1.0]
μ⁰ = [0.0, 1.0, 2.0]

f_sci = SciMLBase.OptimizationFunction(
    (u, p) -> sum(u[i] * (p.μ⁰[i] + log(u[i])) for i in eachindex(u)),
    grad = (g, u, p) -> for i in eachindex(u); g[i] = p.μ⁰[i] + log(u[i]) + 1; end,
    cons = (res, u, p) -> (res .= A * u .- b),
)
u0  = fill(1/ns, ns)
lb  = fill(1e-16, ns)
ub  = fill(Inf, ns)     # SciMLBase requires both bounds as soon as one is given
opt_prob = SciMLBase.OptimizationProblem(f_sci, u0, (μ⁰=μ⁰,);
                                          lb     = lb,
                                          ub     = ub,
                                          lcons  = zeros(1),
                                          ucons  = zeros(1))

sol = SciMLBase.solve(opt_prob, OptimaOptimizer())
println(sol.u)        # ≈ [0.6652, 0.2447, 0.0900]
println(sol.retcode)  # SciMLBase.ReturnCode.Success

[0.6652409557238408, 0.24472847106840737, 0.09003057320775178]
Success

The raw OptimaResult is accessible via sol.original:

raw = sol.original           # OptimaResult
println(raw.converged)       # true
println(raw.iterations)      # number of Newton iterations
println(raw.error_opt)       # final optimality residual

true
28
5.04276355388235e-17

Constraint extraction mechanism

OptimaOptimizer needs the linear constraint matrix and vector explicitly. It extracts them via one of three paths (in order of preference):

1. Explicit parameters: if p is a NamedTuple with fields :A and :b, they are used directly. This is the most efficient path.

2. No constraints: if cons is nothing, the problem is treated as unconstrained (  ,  ).

3. Finite-difference extraction: otherwise, is recovered by forward differencing the constraint function at u0 (one evaluation per species). Accurate for linear constraints; adds extra function evaluations at problem setup time.

To use path 1 directly:

opt_prob_direct = SciMLBase.OptimizationProblem(f_sci, u0, (μ⁰=μ⁰, A=A, b=b);
                                                lb=lb, ub=ub,
                                                lcons=zeros(1), ucons=zeros(1))
sol_direct = SciMLBase.solve(opt_prob_direct, OptimaOptimizer())
println(sol_direct.u)

[0.6652409556346242, 0.24472847109222434, 0.09003057327315153]

Variable scaling details

Internally the SciML interface transforms the problem to scaled coordinates   where   . The scaled conservation matrix   makes the Schur complement well-conditioned across the full range of species concentrations. The transformation is fully transparent: sol.u is always returned in the original unscaled units.

This scaling is particularly important in titration problems where can span 12 orders of magnitude: without scaling, species present at mol effectively vanish from the Newton step and the solver makes no progress toward the correct speciation. ```