Theory

This page describes the mathematical foundations of the Optima solver. The algorithm closely follows the C++ Optima library by Allan Leal (github.com/reaktoro/optima) and the reference by Leal et al. (2014).

Problem statement

We seek the equilibrium composition   minimizing the Gibbs free energy subject to element-mass conservation:

Symbol	Meaning
	mole amount of species
	Gibbs free energy (objective)
	stoichiometric (conservation) matrix
	element abundance vector
	lower bounds (positivity floor  )
	parameter tuple

For an ideal/dilute aqueous solution the gradient and Hessian diagonal are:

Convexity, and what it guarantees

Everything below rests on one structural fact, so it is worth establishing it rather than assuming it.

The ideal mixing term is convex

Write the ideal part of the Gibbs energy, in units, as

on the open positive orthant  . Then

Claim.  , with null space spanned by itself.

Proof. For any  ,

By Cauchy–Schwarz, writing  ,

so   . Equality in Cauchy–Schwarz holds exactly when  , that is  .  

One per phase, and the aqueous convention

A system with several mixing phases contributes one such per phase, each on its own block of variables, so the ideal part of is a direct sum of the Hessians above and is convex whenever each block is. A solid solution is exactly that: the argument does not weaken as phases are added.

The molality convention needs its own completion of the square rather than the same one. With solutes carrying and the solvent , write   and  . Then

the first inequality being Cauchy–Schwarz on the solute block. So the aqueous phase is convex too, and for the same reason.

A pure phase contributes a linear term

A phase of fixed composition has unit activity, so its potential   does not depend on and it contributes to . Hence

which is convex, and is also the source of every numerical difficulty in this solver: the Hessian is singular along those directions by construction, not by accident.

Consequences

With ideal mixing — every activity a mole fraction within its phase, or a molality referred to the solvent — is convex on  , and the feasible set    is a polyhedron, hence convex. Therefore:

1. every local minimum is global, and the minimizer is unique up to the null directions above;

2. the constraints being affine, the linearity constraint qualification holds at every feasible point, so the KKT conditions are necessary and sufficient.

Point 2 is what makes verification possible: a candidate together with multipliers satisfying the KKT system is a proof of global optimality, not evidence for it. Point 1 says that a solver returning different answers from different starting points is not finding different local minima — it is stopping short of stationarity.

Where convexity fails, and what the certificate then proves

Everything above is conditional on ideal mixing, and that condition is not cosmetic.

A non-ideal solid solution carries an excess term. For the Redlich–Kister form used by SolidSolutionPhase,

and a sufficiently positive makes   over an interval. That is not a numerical nuisance, it is a miscibility gap: the non-convexity is the physics, and the phase splits into two of different composition. Aqueous excess terms — Debye–Hückel and the HKF corrections — carry no convexity guarantee either.

Where is not convex, the KKT conditions remain necessary and stop being sufficient. kkt_certificate then proves what it always proves — that the point is a KKT point, feasible, complementary, and with no absent phase supersaturated — and no longer proves that it is the global minimum. Two things follow, and neither is automatic:

• a phase inside a miscibility gap must be allowed to appear twice, as two instances of the same SolutionPhase with different compositions, or the model cannot represent the answer at all;

• the tangent-plane test admits a phase from outside; it does not test whether an already-present phase would rather split. Deciding that is Michelsen's stability problem on the phase itself, and it is not what this solver does.

Every mixing model this solver has been exercised on is ideal, so the certificate has meant global optimality throughout — but that is a property of those models, not of the method, and it should be checked before it is relied on.

Log-barrier interior-point method

Barrier augmentation

To enforce   strictly without explicit inequality constraints, a logarithmic barrier term is added to the objective:

The barrier weight   is driven to zero over an outer loop; as   the barrier minimizer converges to the original constrained minimizer.

KKT conditions

At a stationary point of the barrier-augmented Lagrangian

the first-order KKT conditions are:

where division by   is component-wise. KKTResidual stores together with and .

The optimality error, and why the obvious one cannot work

Let     denote the slacks and

the gradient of the Lagrangian without the bound term. The barrier stationarity condition is

The residual form diverges at the bounds

Measuring   as it stands, i.e. by   , is unusable, and not marginally so.

Claim. Along any sequence approaching a solution at which species is absent ( ), is unbounded for a fixed relative error in .

Proof. Let    with the exact slack and the relative error. At the exact point  , so

for any fixed  .  

Numerically: with   and   — an ordinary state of affairs for a mineral that is not present — the term alone is . On a cement equilibrium this quantity started at   and never fell below  , so tol = 1e-4 failed exactly as tol = 1e-10 did: no tolerance was attainable, and the failure was not one of accuracy but of the criterion itself.

Two things followed, both silent. The barrier update should_reduce_barrier requires the inner loop to reach before may fall, so stayed at its initial for the whole solve; and every solve exhausted max_iter and returned whatever iterate it had reached, while the feasibility error   was meanwhile reaching .

The guard intended to prevent this — excluding variables within   of their bound — could not fire either: with   that threshold is , so nothing was ever excluded.

The complementarity form is bounded and equivalent

Introduce the bound multipliers   explicitly. The KKT system of the barrier subproblem is

Dual feasibility defines  ; substituting it into complementarity gives the equivalent scalar conditions

Claim.   and    have the same solutions, and   is bounded where is not.

Proof. Since   strictly at any interior point,   vanishes iff does, which gives the equivalence. For the boundedness, is finite at any interior point and is bounded above by the total amount, so    ; whereas the computation above shows    as   at fixed .  

This is the measure Ipopt reports, with the bound multiplier eliminated by dual feasibility rather than carried as a variable (Wächter & Biegler 2006, §3.5). On the same cement equilibrium it reads   where the residual form read   — fourteen orders of magnitude, on the same iterate.

The feasibility error is not rescaled:   is in moles and already means something absolute.

Convergence criterion

What this method does not do

Making the error measure meaningful does not make the iteration converge, and it is honest to separate the two.

With the error now readable, the trace on a cement equilibrium shows the step capped by the fraction-to-boundary rule of Fraction-to-boundary step limit at every iteration — equal to throughout, falling below — so convergence is linear and tol is not reached. The mechanism is the one identified above: a pure phase has   exactly, so its Newton direction is governed by the barrier term alone and is large, and the boundary cuts it.

Three remedies were implemented and measured, and all three are rejected because they moved the chemistry the wrong way:

remedy	effect on the KKT error	effect on the answer
reduce on stalling	barrier finally falls	worst element imbalance   mol
project each iterate onto	feasibility	pore solution pH
primal-dual with explicit	feasibility to	aluminate assemblage lost entirely

The structural fix is not a modification of this iteration but a different formulation: solving the KKT system directly in the space of element potentials, where the aqueous species are parameterized by so that positivity is automatic and the fraction-to-boundary rule has nothing to act on. That is implemented in ChemistryLab.DualEquilibriumSolver, which uses this solver to reach a neighborhood and then certifies the result.

Newton on the KKT system in multiplier space

The interior-point method above minimizes f by walking the interior. This section describes the alternative the package also provides, dual_newton_solve, which solves the KKT conditions directly. It is the Brinkley–Karpov formulation of the geochemical Gibbs-minimization codes, stated here for the general convex program.

The formulation

Let u := -A^\top y. The stationarity conditions split by the nature of the variable:

The second line is a stability criterion: a bounded variable is positive exactly when u_i - g_i vanishes, and zero when that index is negative. In a chemical system it is the statement that a phase is present iff it is saturated.

Why it is well conditioned

Parameterizing the interior variables by   makes their positivity automatic. The fraction-to-boundary rule of Fraction-to-boundary step limit — which caps the interior-point step at every iteration on a cement equilibrium — therefore has nothing to act on for them. Only the bounded variables carry a bound, and they are handled by an active set, exactly.

This is not the log reparameterization that variable_space = Val(:log) performs, and the distinction matters. Composing the objective with exp gives

which is negative wherever    — for a chemical potential of order , everywhere. f ∘ exp is not convex. Here the logarithm is applied to the KKT equations, solved as a square nonlinear system; convexity of the original problem is what makes that system's solution unique.

The reference variable

One interior variable may have h_i bounded above, so that its condition h_i = u_i - g_i has no solution for an arbitrary y. In a chemical system that is the solvent: its activity is a mole fraction, so   always. Inverting it is not merely slow, it can be infeasible, and an inner loop that included it can never report convergence — which then invalidates the outer Jacobian, since that Jacobian is derived on the assumption that the inner conditions hold exactly.

Such a variable is declared through j_ref and carried by the outer system, where the equality constraints determine it.

Phases with no solvent

The paragraph above is written for a phase that has a solvent. A solid solution has none: every member is a mole fraction, so every h_i is bounded above by zero, and the argument that singled out one variable now applies to all of them at once. Inverting h_i = u_i - g_i member by member is then impossible whenever the right-hand side is positive, and an iteration that tries answers by sending the member's amount to infinity.

Fixing the reference's amount and inverting the others does not repair it. With x_\text{ref} held and ideal mixing, the fixed point gives    with   , hence

which has no positive solution once   — precisely the situation in which the phase is supersaturated at the current y and the reference cannot hold it back. Measured on a cement, the four C-S-H end-members all ran to the clamp at    mol, and the outer Jacobian was then computed on that.

What the potentials do determine, always, is the composition. Stationarity of every member gives  , so

with the phase total N as the outer unknown. Both are evaluated with the maximum factored out, so nothing overflows for any y. The phase equation is the same expression as the tangent-plane measure below, which is what one would want: a phase is admitted and held stationary by one quantity, not by two that could disagree.

SolutionPhase carries this as mole_fraction. It must stay false for an aqueous phase: there the solutes' activities are molalities, unbounded above, and their dependence on u is exactly what determines the multipliers — remove it and y is barely determined at all.

Degenerate components

A row k with   need not be degenerate. With  ,

only if the non-zero entries of the row share a sign — a sum of non-negative terms vanishes only term by term. A row with entries of both signs permits cancellation and forces nothing. degenerate_components implements exactly this test.

The distinction is not academic: the H+ row of a chemical system carries +1 for H+ and −1 for OH-, so its zero total is the ordinary state of pure water. Declaring it degenerate removes the entire acid–base system and returns pH 7.000 with the solid undissolved.

Where a row is degenerate its multiplier is determined by nothing and the Jacobian is singular in that direction. The row is then replaced by  , which keeps the system square and makes the pinning of those variables consistent with their stationarity rather than merely imposed on it.

Termination

The outer loop admits one variable per round, the most violated, and records the active sets it has visited. Since there are finitely many subsets and each round either terminates or visits a new one, the loop terminates. Admitting a batch instead feeds a cycle in which a variable is admitted, driven negative, dropped, and readmitted — observed on a cement without limestone as a solve converged to 2e-12 of the wrong subproblem, with an excluded variable violated by 10.9.

The starting multipliers

A least-squares fit of     minimizes a stationarity residual and says nothing about b, so the potentials it produces are consistent with no composition in particular. The last candidate is therefore solved for instead of fitted. Treating every species as an ideal one whose activity is its own amount, the Lagrangian minimizes in closed form,  , and the dual becomes

is smooth and strictly concave, so Newton with a backtracking line search converges from any starting point — there is no active set to guess and nothing to stall in. This is Brinkley's method (White, Johnson & Dantzig, 1958), and its solution is the y for which the ideal composition conserves matter exactly.

It is tried last, not first. A caller replaying a trajectory hands over a composition that is nearly the answer, the fits built from it converge on the first attempt, and this solve is then never run; ordered the other way it costs a warm-started replay three digits of element balance for no gain. The attempts are ranked by the KKT error each reaches, so the answer returned is the best of them whether or not any crossed the tolerance.

What an active set must satisfy, and how it is searched

Two conditions decide whether an active set can support a solution at all, and both are consequences of the KKT system rather than choices.

A bound-constrained variable held active imposes u_i = g_i, i.e.    — one linear equation in y. A mole-fraction mixing phase imposes  , one more, nonlinear but still a condition on y alone. A phase with a solvent does not, since its reference equation involves the composition. Hence

which is Gibbs' phase rule in the form the dual system takes, and in addition    — two dependent columns demand a fixed relation between their g_i that no database satisfies, which is exactly what two polymorphs of one composition amount to. An active set breaking either condition has no solution: the residual cannot reach zero for any iterate, and the least-squares step spreads the violation over the rows instead of removing it.

Complementarity supplies the leaving rule. For a bound-constrained variable the conditions are  ,  ,   with    . Testing only   covers one half: a variable held active while undersaturated violates complementarity as plainly, and no Newton iteration repairs it, because its own equation   is the one that cannot hold. That test is meaningful only on a point that solves the current subproblem — while the inner Newton is still working, is a transient.

The search over sets is a descent method, and the quantity it descends is the KKT error of the whole problem — stationarity, element balance, and the worst violation among the phases held absent — not the residual of the subproblem the current set defines. The distinction decides the outcome: a set that omits a phase the solution needs solves its own equations exactly, so ranking states on that residual rewards leaving phases out. The best state visited is what the solver returns.

Exchange, not rejection

An admission that fails to converge is not evidence against the variable admitted. The inner Newton stops the moment an active variable falls below its bound, and that is the departing variable announcing itself: the entrant stays, the drop list removes the other, and the two together are the active-set exchange.

Reading that failure as the entrant's fault is what an earlier version did, and the consequence was silent. On a cement, ettringite and monosulphate compete for the same sulfate, so admitting one necessarily drives the other out; rejecting ettringite permanently let the solve converge — to 2e-12 stationarity and 5e-12 element balance — onto an assemblage in which ettringite was absent and supersaturated by 14.8. Only the certificate caught it.

The entrant is reconsidered only when the Newton failed with nothing leaving, which is the genuine over-determination this guard was written for: two bound variables declared stationary whose formulas are dependent modulo the mixing phases, where the residual cannot reach zero at all. Even then the veto lasts only as long as the context that produced it, and a vetoed variable that is still supersaturated prevents the run from being reported as converged — the candidate list is filtered, the KKT conditions are not.

The certificate

kkt_certificate checks the conditions at any point, whatever produced it. For a convex program they are sufficient, so a certificate is a proof of global optimality.

Two splits decide whether the check means anything. A variable at its bound obeys the inequality, not the equality: imposing the equality on an amount held at 1e-16 whose stationarity value is e^{-300} misstates h_i by 263 units, and the check then reports a residual of 74 for a point solved to 5e-12. And a variable carrying a degenerate component is excluded from both tests, for the reason above.

Newton step via Schur complement

KKT linear system

At each Newton iteration we solve the     saddle-point system:

where   is the barrier-augmented Hessian diagonal:

Schur complement reduction

Because is diagonal, can be eliminated analytically. From the first block row:    . Substituting into     gives the   Schur system:

Once is found, is recovered by back-substitution.

Implementation. is built as a single BLAS GEMM:   where   is computed in-place. The RHS is then the BLAS GEMV   , and is recovered by the BLAS GEMV . All three operations reuse the same pre-allocated buffer (field AoverH of NewtonStep).

The total cost is to build and to factor it. Since (number of conserved elements) is typically , this is far cheaper than factoring the full  -dimensional system.

Numerical conditioning

Two techniques stabilize the Schur solve when some conservation rows correspond to absent species (zero element budget):

1. Tikhonov regularization: add with   before factorization, preventing near-zero pivots.

2. Diagonal equilibration: scale row and column by so all diagonal entries equal 1, reducing the condition number from to in titration-type problems.

Conservation matrix canonicalization

Canonicalizer decomposes via QR with column pivoting:

The LU factorization of is cached and reused across Newton steps, reducing each back-substitution to rather than . When is fixed across a sequence of solves (e.g. a temperature scan), pass the pre-built Canonicalizer to solve to skip the QR entirely.

The rank and the basis order are two questions, and each needs its own answer.

For the rank, the column scaling is noise:    for any positive s. But the pivoted-QR test compares each pivot to the largest, and the SciML interface scales columns by each variable's starting value; warm-starting from a converged equilibrium spreads those over ten orders of magnitude and the smallest honest pivot falls below the tolerance. The rank came out one short, B was built with m-1 columns, and the failure surfaced as a "matrix is not square" error from LAPACK. It is therefore read off a column-equilibrated copy.

For the basis order, that same scaling is exactly the information wanted. The null-space step asks the basic variables to absorb the infeasibility, through   , so they must be the ones that can move — the abundant species, not a trace ion pinned at its bound. Equilibrating before pivoting throws that away: on an LC³ equilibrium the basis then held species at 1e-16, the particular solution asked them for 1e4 mol, and the dual step came back at 1e31. The order is taken from the matrix as given, which is how Optima prioritizes its own basis.

A genuinely rank-deficient A — a conservation law that is a combination of the others — is reported as such instead of being discovered inside a factorization.

The reduced Hessian must be equilibrated

With the canonical partition, the step may be taken in the null space of A:

h is the barrier-augmented curvature  . In a chemical system the amounts span ten orders of magnitude, so h spans twenty and more — on a cement it ran from 2.5 on the solvent to on a species at its bound. The reduced Hessian inherits that spread, its condition number passes anything Float64 can carry, and the Cholesky then succeeds while returning a direction that is noise:   ,   , and NaN two iterations later.

Scaling by the square root of the diagonal is exact — with  , solving    returns the same — and it normalizes every diagonal entry to one. The Schur-complement branch already did this; the null-space branch, which is the default, did not.

Convergence, and the level it is measured at

The barrier subproblem's stationarity is sᵢ (∇f + Aᵀy)ᵢ = μ, so

vanishes at the solution of that subproblem, for any μ. It is the right quantity to decide when to tighten the barrier, and the wrong one to decide that the problem is solved: a point satisfying it at μ = 10^{-4} is O(10^{-4}) from the optimum. Convergence is therefore declared on

the same residual at μ = 0 — Ipopt's E_0, Wächter & Biegler (2006), Algorithm A, step 2. It is the true KKT error and cannot be satisfied at a loose barrier.

That test only becomes attainable with a barrier schedule the inner Newton can keep up with. At κ_μ = 0.1 the barrier outruns it and E_0 plateaus just above the tolerance without crossing it — 312 iterations to reach 9.999\times10^{-11} against 10^{-10}. At Ipopt's κ_μ = 0.2 the same problem converges in 30.

Feasibility of the starting point

The barrier method must be started on    , not near it, and the reason is the line search. The filter is bypassed only when the current point is feasible; while it is not, a step is accepted either because the constraint violation drops by a relative — unreachable once the fraction-to-boundary limit is itself below — or by Armijo along a direction partly spent restoring feasibility, which need not be a descent direction at all. Measured on an LC³ equilibrium, a start carrying     had all forty trial steps refused at every barrier level, and the solve reported MaxIters on the point it began at.

Positivity is enforced first and feasibility after, so the clamp cannot undo the projection. Two exact routes are then tried in order:

1. Solve for the basic amounts given the others,   . One triangular solve on a factorization that already exists, and this is how Optima does it.

2. Non-negative least squares on the slacks   , by Lawson–Hanson. Needed because a component total can be negative — the H⁺ row of a cement is mol — and no basis of abundant species can produce it. NNLS terminates finitely and its residual is zero whenever the budget is attainable, which it is whenever the budget came from a real composition.

NNLS returns a solution supported on at most variables, the rest at exactly their bound with zero slack — unusable as a barrier start, since the first negative step component gives  . They are lifted to the slack the barrier itself would give them: a variable held at its bound settles at  , about for the initial   with chemical potentials of order – in RT units. The matter that adds is then removed exactly, on the support, by a minimum-norm correction, so the point is feasible to machine precision and strictly interior.

Alternating projections onto the two sets converge to the same place in theory and far too slowly to be usable — on that same problem they reached   and then advanced by less than a tenth of a percent per sweep. They remain only as a last fallback.

Fraction-to-boundary step limit

Before the line search, the full Newton step is scaled to keep all components strictly above their lower bounds:

Additionally, unstable variables — species that are near their lower bound (   ) with   (gradient pushing toward the bound) — receive a further reduced step:

This prevents numerical oscillations when a species is in the process of precipitating or dissolving completely.

Filter line search

The line search follows Wächter & Biegler (2006). A filter is a Pareto set of pairs    ; a new point is acceptable to the filter if it is not dominated by any entry already in the filter.

Starting from  , the algorithm backtracks with factor   until the candidate       satisfies:

• Filter acceptance: not dominated by any entry in the current filter, and

• Sufficient decrease on the barrier objective (Armijo condition):

or a sufficient feasibility decrease:

When the current iterate is already feasible ( ), the filter is bypassed and only the Armijo condition on is checked, switching the method to a pure descent algorithm for the final convergence phase.

Outer barrier loop

The outer loop reduces on a geometric schedule:

The inner Newton loop for each fixed runs until the KKT error satisfies

so the inner tolerance tightens automatically as  , avoiding unnecessary Newton iterations in the early (exploratory) phase.

Sensitivity analysis

At convergence , the implicit function theorem applied to the KKT system    gives

Two parameter families are of direct chemical interest:

Response to element budgets

The right-hand side for perturbation of is    , giving:

Sanity check: summing over species ,   — the extra mole of element is fully redistributed among the species.

Response to standard potentials

The right-hand side for perturbation of is   , giving:

Implementation. Both sensitivity matrices are computed with a single batched solve (BLAS TRSM) followed by a BLAS GEMM, rather than sequential scalar solves:

where   (the same buffer built during the last Newton step). Using a matrix right-hand side triggers BLAS level-3 (TRSM + GEMM) instead of level-2 (TRSV + GEMV) calls — a significant speedup when  .

The total cost is   — negligible compared to the solve itself.

The warm-start cache, and the caller's initial point

OptimaOptimizer carries a cache of its previous solution, and with warm_start = true it may reuse it. The semantics matter enough to state precisely, because the natural implementation is wrong.

Let the caller supply an initial point . The cache holds , the solution of whatever problem this algorithm object last solved — not necessarily this one, since the object is reusable and reused.

The rule. The cache is consulted only when

1. it has the same dimension as the current problem, and

2. carries no interior information, i.e. every component sits at its lower bound,    for all .

Otherwise is used as given.

Condition 2 is the substantive one. Starting from regardless discards an initial point the caller chose deliberately, and the consequence is not academic: a chemical-kinetics run re-speciating at every accepted step leaves its final composition in the cache, so replaying the same trajectory through the same algorithm object started every solve from the 28-day state. On an ordinary Portland cement that returned a pore solution at pH 14.2 with 0.31 mol of ettringite and no monosulphate, where honoring the caller's guess gives pH 12.58 with the sulfate entirely in monosulphate — the same trajectory, the same constraints, the same guess. It also silently defeated the caller's own warm-start logic during the run.

The cache is therefore a convenience for repeated solves where the caller has nothing better to offer, and never an override.

Variable scaling in the SciML interface

The OptimaOptimizer SciML interface automatically scales each species by its starting value  :

This transforms the scaled problem so all   at the starting point, making the Schur complement well-conditioned across the multi-decade concentration ranges typical in chemical speciation (e.g. pH 1–13 where varies over 12 orders of magnitude). The scaling is transparent: the returned solution is always in the original units.

References

• Allan Leal, Optima — C++ library for chemical equilibrium optimization, ETH Zürich. github.com/reaktoro/optima

• Leal, A.M.M., Blunt, M.J., LaForce, T.C. (2014). Efficient chemical equilibrium calculations for geochemical speciation and reactive transport modelling. Geochimica et Cosmochimica Acta, 131, 301–322. https://doi.org/10.1016/j.gca.2014.01.038

• Wächter, A., Biegler, L.T. (2006). On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathematical Programming, 106(1), 25–57. https://doi.org/10.1007/s10107-004-0559-y