A geometric reformulation of the bilevel parameter optimization problem to a single level non-linear programming problem with applications to phase equilibria

Abstract: Phase equilibrium problems are central to chemical engineering, underpinning tasks ranging from separation process design to the development of thermodynamic models. A particularly challenging computational task is the generation of phase envelopes: rigorously fitting thermodynamic models to real world data with well behaved predictions requires solving a computationally expensive bilevel optimization problem. We present a geometric reformulation of this parameter estimation problem that restates the bilevel program as a single level problem that is significantly easier to solve. The solution of the single level problem is proven to be the globally optimal solution of the bilevel problem when specialized global optimization solvers are used. In addition, the method retains the constraints that guarantee a well behaved fitted model, such as enforcing the correct number of phase splits, excluding spurious phases, and ensuring stability in regions of instability. This allows the practitioner to reliably and efficiently fit mathematically complex thermodynamic models to data, and potentially enables highly accurate and rigorous modelling of problems in computational thermodynamics that were previously intractable. Finally, an algorithm is presented that is proven to converge for any black box thermodynamic model. Only an expression of the Gibbs free energy is required, no derivatives are needed, and convergence is guaranteed for the broadest class of non-smooth, non-continuous models.

• The paper replaces nested Gibbs-energy minimization and stability checks with a dual-manifold non-negativity condition and a single-level SHGO-based nonlinear program, with a zero objective certifying globally consistent phase-equilibrium fits.
• The method matches or improves published binary NRTL benchmarks, reproduces a CEOS result within 1.65% while using roughly two orders of magnitude less reported computation time, and identifies localized defects when models are misspecified.
• Industrial regressions predicted held-out bubble pressures with average absolute deviations of 9.74% for sour gas, 1.97% for carbon capture, and 1.27% for R-407C, while exposing limitations near critical regions and for weakly identified parameters.

The paper develops a geometric reformulation of the bilevel parameter-estimation problem that arises when thermodynamic models are fitted to phase-equilibrium data, and demonstrates it on binary liquid–liquid equilibrium (LLE) benchmarks and three industrial multicomponent case studies. In the classical formulation of Mitsos, Bollas and Barton, the upper level minimizes a model–experiment mismatch over model parameters $\mathbf{p}$ while the lower level requires, for every candidate parameter set, a global Gibbs free-energy minimization plus a tangent-plane stability check. That structure is computationally expensive: deterministic solutions via branch-and-bound over KKT-reformulated single-level programs, as in the cubic-equation-of-state (CEOS) extension of Glass, Djelassi and Mitsos, cost roughly 4.5 hours per binary case study, a cost the original authors themselves identified as prohibitive for industrial-scale problems. The present work replaces the nested lower-level programs with a single geometric object — the dual manifold — and solves the resulting single-level nonlinear program with the derivative-free, deterministic simplicial-homology global optimization (SHGO) algorithm.

The dual manifold construction

The central object is defined on the composition simplex $\Delta^n$. For each experimental tie-line $i$ with endpoint compositions $\mathbf{x}^{\alpha,i}$ and $\mathbf{x}^{\beta,i}$, an affine chord $h_i$ interpolates the reduced Gibbs energy $g = G/(RT)$ between the two anchors. The composite supporting hyperplane is the pointwise minimum $h(\mathbf{x}) = \min_i h_i(\mathbf{x})$, and the dual manifold is the shifted surface

$$\mathcal{M}(\mathbf{x};\mathbf{p}) = g(\mathbf{x},\mathbf{p}) - h(\mathbf{x}).$$

By construction $\mathcal{M}$ vanishes at every tie-line endpoint. A parameter vector reproduces every datum as a globally stable equilibrium under the Gibbs–Michelsen tangent-plane criterion if and only if $\Delta^n$0 over the entire simplex; a negative excursion is a topological defect — a parameter set that interpolates the data points yet violates global stability elsewhere, precisely the mechanism behind documented spurious-phase failures of standard regressions such as the CF$\Delta^n$1/CHF$\Delta^n$2 SRK fit and the acetone/water Peng–Robinson fit.

The single-level objective $\Delta^n$3 aggregates five non-negative components: a tangent-plane sub-sampling residual evaluated at four bracketing points per tie-line, a plane-violation term over the set of defective minima $\Delta^n$4 located by an inner SHGO call on $\Delta^n$5, a topological pull of each defect toward its nearest equilibrium endpoint, a local-topology amplitude measuring the joint composition–energy extent of each negative connected component, and a phase-assignment partition relevant to multi-root equation-of-state models. A two-stage solver runs the cheap tangent-plane pass first (no inner global search), invokes the defect-aware stage only when the inner SHGO call detects surviving defective minima, and loops to a fixed iteration limit.

Theoretical guarantees

Three proof obligations are discharged in the appendices. First, the equivalence $\Delta^n$6 on $\Delta^n$7 if and only if $\Delta^n$8 is a globally consistent fit is proved for continuous $\Delta^n$9 under an explicitly stated data-consistency assumption on the tie lines, combining the Gibbs–Michelsen tangent-plane criterion with the Mitsos–Barton dual extremum principle. Second, a compactness argument establishes that $i$0 is a sufficient condition for global stability at every datum, independent of the inner solver's behavior. Third, convergence of the two-stage outer loop is split into an unconditional and a conditional half. Whenever a parameter vector with $i$1 is returned, its global optimality and thermodynamic consistency are certified a posteriori by the non-negativity of $i$2 alone. That such a point is found whenever it exists is inherited from the deterministic adequate-sampling theorems of SHGO, which are finite and non-probabilistic. The one step that remains an assumption is that the bifurcation strata of $i$3 carry zero Lebesgue measure; this is justified for real-analytic model families (NRTL, Wilson, UNIQUAC, cubic equations of state) but asserted rather than proved for genuinely tabulated black-box $i$4. Notably, the construction is model-agnostic: only an expression of $i$5 is required, no derivatives, no KKT system, and no cubic equality constraint — root selection for CEOS models is replaced by the convex envelope of the $i$6-branches, which removes the root-discrimination prerequisite that Glass et al. identified as the obstacle to extending their formulation to non-cubic equations of state such as PC-SAFT.

Validation on the LLE benchmark

The four binary NRTL LLE systems of Mitsos, Bollas and Barton serve as the direct literature comparator, with error metrics matching their published table. On Case 1 (n-butyl acetate–water) the bilevel fit attains an absolute residual of $i$7 against the published $i$8, with per-phase relative errors two to three orders of magnitude smaller. Case 3 (n-octanol–water) achieves $i$9, two orders of magnitude below the literature value, with full 50/50 phase-envelope coverage. Cases 2 (n-butanol–water) and 4 (furfural–2,2,5-trimethyl-hexane) show larger absolute residuals ($\mathbf{x}^{\alpha,i}$0 and $\mathbf{x}^{\alpha,i}$1); the authors trace these to the common-tangent phase-envelope evaluator losing its convergence basin near the consolute boundary — only 35/50 and 46/50 grid temperatures return a two-phase tangent even at the published literature parameters — rather than to the bilevel procedure. The per-temperature dual-manifold record is reported candidly: Cases 1 and 3 are defect-free, while Cases 2 and 4 retain negative lobes (up to $\mathbf{x}^{\alpha,i}$2 at 393.15 K for Case 2), which the authors present as evidence that binary NRTL is misspecified for those mixtures, not as defect-free fits. The fitted temperature-independent coefficients $\mathbf{x}^{\alpha,i}$3 agree with the literature to within a few percent in most cases, while the $\mathbf{x}^{\alpha,i}$4 and $\mathbf{x}^{\alpha,i}$5 coefficients diverge sharply, a known degeneracy of NRTL temperature dependence over narrow temperature ranges.

Industrial case studies

Three multicomponent systems are regressed on binary subsystems only, with the multicomponent behavior predicted by extrapolation and every tie line certified against the dual-manifold criterion. For the sour-gas ternary CH$\mathbf{x}^{\alpha,i}$6 + CO$\mathbf{x}^{\alpha,i}$7 + H$\mathbf{x}^{\alpha,i}$8S with Peng–Robinson, the fitted CO$\mathbf{x}^{\alpha,i}$9/H$\mathbf{x}^{\beta,i}$0S interaction parameter reproduces the industrial reference to within 0.05%, and the model predicts 31 held-out ternary bubble pressures at AAD($\mathbf{x}^{\beta,i}$1) = 9.74%; 27 of 31 tie-line sections are defect-free, the four exceptions being near-critical sections with shallow lobes of magnitude at most $\mathbf{x}^{\beta,i}$2, which forcing the defect stage to suppress would cost 12.5% in bubble-pressure deviation — a diagnostic of structural Peng–Robinson misspecification near the critical locus rather than of the regression. For the carbon-capture quaternary CO$\mathbf{x}^{\beta,i}$3 + N$\mathbf{x}^{\beta,i}$4 + Ar + O$\mathbf{x}^{\beta,i}$5, out-of-sample ternary validation yields AAD($\mathbf{x}^{\beta,i}$6) = 1.97% over 18 held-out bubble pressures with all 18 sections defect-free; the fitted CO$\mathbf{x}^{\beta,i}$7/Ar and CO$\mathbf{x}^{\beta,i}$8/O$\mathbf{x}^{\beta,i}$9 parameters ($h_i$0, $h_i$1) fill a genuine gap, since reference values were assumed zero. For the zeotropic refrigerant blend R-407C, the 44 held-out ternary bubble pressures are reproduced at AAD($h_i$2) = 1.27% with an unqualified stability record: all sections satisfy $h_i$3 with the worst interior minimum at machine precision ($h_i$4), demonstrating that the certificate $h_i$5 is genuinely attainable when the model is adequate. The authors note honestly that for this near-ideal blend the objective is flat in the interaction parameters, so the fitted values are weakly identified.

A separate appendix replicates the Glass et al. CEOS proof-of-concept on the C$h_i$6H$h_i$7/H$h_i$8S system, recovering their published interaction parameters to within 1.65% (Peng–Robinson) and 0.55% (SRK), and reproducing both of their published figures.

Computational cost and limitations

The cost comparison is stated with appropriate caveats. The CEOS replication converges in 164.7 s (PR) and 106.9 s (SRK) of wall clock on a single commodity core against the approximately 4.5 hours per case study reported by Glass et al. with BARON — a factor of roughly one hundred, which the authors explicitly decline to inflate: hardware differs, CPU time and wall clock are different quantities, and a branch-and-bound $h_i$9-optimality certificate is not the same work as a sampling-based inner call. The four LLE fits take between five and twenty minutes each. The structural asymmetry that survives the caveats is that the method never encodes the cubic equality as a constraint and never forms the lower-bounding subproblems identified as the bottleneck in the branch-and-bound approach.

Several limitations are conceded at the point where they bear on the results. The bifurcation-strata assumption is unproven for arbitrary black-box $g = G/(RT)$0, and a proof for tabulated $g = G/(RT)$1 is left open. Adequate sampling of the inner call is certified only when the returned minimiser pool meets the Gibbs phase-rule bound on the number of coexisting phases; where it does not, the reported defect counts are lower bounds on defects actually present, which is how the supplementary stability archives should be read. In every industrial case the defect-aware stage acted as a certifier rather than a refitter: it located, counted and reported defective sections but did not drive the parameters to $g = G/(RT)$2 where the underlying model is inadequate, so the certificate is informative in both directions — unconditional when attained, and a quantified, localized failure report when not. The common-tangent envelope evaluator, not the bilevel method, limits Cases 2 and 4, and a more robust construction (homotopy continuation or Lagrangian duality) is deferred to future work. No fitting result is claimed for the software's remaining model adapters (Wilson, UNIQUAC, UNIFAC, SRK, PR78, TWU, PC-SAFT), which are capabilities of the released code rather than evidence in the paper; carrying the model-agnostic argument through to an actual PC-SAFT regression remains an open question, as does certification that no NRTL or UNIQUAC parameter set can reproduce the island-type ternary LLE topologies analyzed by Olaya et al., which the dual manifold makes tractable but which is not demonstrated here.

Conclusion

The paper recasts the bilevel thermodynamic parameter-estimation problem as a single-level nonlinear program over a dual-manifold non-negativity condition, with proofs that a vanishing objective certifies global thermodynamic consistency and that the two-stage SHGO solver terminates finitely under stated sampling and regularity assumptions. Empirically, the method matches or improves published LLE benchmark fits, replicates an established CEOS bilevel result in about two orders of magnitude less time, and delivers binary-only regressions with certified stability extrapolated to three industrial multicomponent systems at 1.27–9.74% bubble-pressure error. The certificate it issues is strongest when the model is adequate and degrades gracefully into a localized defect report when it is not.