the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 04 Nov 2024
An unconstrained formulation for complex solution phase minimization
Abstract. Prediction of mineral phase assemblages is essential to better understand the dynamics of the solid Earth, such as metamorphic processes, magmatism and the formation of mineral ore deposits. While recently developed thermodynamic databases allow the prediction of stable phase mineral assemblages for an increasing range of pressure, temperature and compositional spaces, the increasing complexity of these databases results in a significant increase of computational cost, hindering our ability to perform realistic models of reactive fluid/magma transport. Presently, prediction of stable phase equilibrium in complex systems is therefore largely limited by how efficiently single phase minimization can be performed, as more than 75 % of the total computational time is generally dedicated to individual solution phase minimization. This limitation becomes critical for non-ideal solution phase models that involve both a large number of chemical components, and mixing on a large number of sites, resulting in many inequality constraints of the form 0 ≤ xMl ≤ 1, where xMl is the fraction of element l mixing on site M.
Here, we present a general reformulation of complex non-ideal solution phases from the thermodynamic database of Holland et al. (2018), which comprises equations of state for multiple mineral solid solutions appearing in magmatic systems, as well as multicomponent silicate melt and aqueous fluid phases. Using a nullspace approach, inequality constraints governing the site fractions are transformed into equality constraints, and the resulting problem is turned into an unconstrained optimization problem, subsequently optimized using efficient gradient-based methods. To test our formulation, we apply it to several equations of state for solution phases known for their complexity and compare the results of our approach against classical optimization algorithms supporting inequality constraints.
We find that the BFGS algorithm yields by far the best performance and stability with respect to the other investigated methods, improving the minimization time of individual solution phase by a factor ≥ 10. We estimate that our new approach can improve the computational time of stable phase equilibrium by a factor ≥ 5, thus potentially allowing to model realistic reactive fluid/magmatic systems by directly integrating phase equilibrium calculations in multiphase thermomechanical codes.
Competing interests: Boris Kaus is editor of the special issue
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this preprint. The responsibility to include appropriate place names lies with the authors.- Preprint
(2960 KB) - Metadata XML
- BibTeX
- EndNote
Status: open (extended)
-
RC1: 'Comment on gmd-2024-197', Anonymous Referee #1, 15 Jan 2025
reply
The authors present a modification of their code MAGEMIN to allow for more efficient Gibbs free energy minimization. They state that the advantages of this modification arise from the new minimization algorithm being “unconstrained” as opposed to the “constrained” algorithm used previously. This distinction is misleading and the authors should change the paper to more clearly state the change in strategy. Simply put, the new algorithm is NOT “unconstrained”.
The central problem is simply stated: minimize the Gibbs free energy of a phase as a function of the fractions of cations X_i, I=1,n on one or more crystallographic sites. The minimization is subject to two constraints:
1) The cation fractions must sum to 1:
sum_i^n X_i = 1
and
2) they must be non-negative
X_i >= 0
Constraint (1) can be implemented either by performing a constrained minimization for which (1) forms an auxiliary statement to the minimization problem OR by minimizing over the null space of the constraint (1). The latter approach is appealing because it reduces the dimensionality of the problem (by the number of crystallographic sites), and avoids the auxiliary statement of constraint.
The primary focus of this paper is on the implementation of the null space approach. And it is on the basis of the removal of the auxiliary statement (1) that they describe their modified algorithm as “unconstrained”.
However, it is NOT unconstrained. The reason is constraint (2). This must still be dealt with: the minimization must still be subject to the bounds described by constraint 2. Thus the minimization is still constrained, not unconstrained. In other words, they have removed an equality constraint, but NOT the inequality constraint.
Of course the authors DO end up applying the inequality constraint, although this application is somewhat buried in the details (Eq. 32).
My argument is not with the method itself, nor primarily with the claims of greater efficiency (although see below). It is with what I think is a misleading description of the modification. I implore them to characterize their modification more clearly and correctly.
Some further comments.
1) Readers may get the mistaken impression that the null space approach is new in the context of petrological Gibbs free energy minimization codes. It is not, and an example of another minimization code that uses the null space approach is HeFESTo (Stixrude and Lithgow-Bertelloni, 2011). Nor do the authors of HeFESTo claim priority for the null space approach, but it is the closest implementation known to this reviewer in terms of intended application (i.e. to petrology).
2) The claimed superiority of their new implementation to SLSQP is puzzling to this reviewer. One reason for the source of puzzlement is that HeFESTo uses SLSQP and finds a very high (essentially perfect) rate of success and precision as documented in Stixrude and Lithgow-Bertelloni (2021). Perhaps the reason is that HeFESTo uses the null space approach for the equality constraints and relies on the constraint facility of SLSQP only for the inequality constraint. Maybe in the current paper, the authors are instead relying on SLSQP to take care of the equality AND inequality constraints. Therefore, to avoid confusion in the literature, I suggest the following test:
Perform the SLSQP minimization(s) again, but by combining SLSQP with the null space approach.
3) A final comment on SLSQP. The authors state that SLSQP sometimes fails because it violates inequality constraints and then cannot return to the feasible space. But this should not be true. According to NLOPT documentation SLSQP is guaranteed to respect inequality constraints at all intermediate steps of the minimization. Perhaps the problem is that SLSQP sometimes does venture into the space where one or more X_i are exactly zero. If this is the case, it is a problem easily solved, by setting the inequality constraint instead to:
X_i >= epsilon
very much like their own implementation of the inequality constraint.
Citation: https://doi.org/10.5194/gmd-2024-197-RC1 -
AC1: 'Reply on RC1', Nicolas Riel, 27 Jan 2025
reply
We thank the reviewer for reading our manuscript and giving helpful comments that will allow us to clarify our work.
The authors present a modification of their code MAGEMIN to allow for more efficient Gibbs free energy minimization. They state that the advantages of this modification arise from the new minimization algorithm being “unconstrained” as opposed to the “constrained” algorithm used previously. This distinction is misleading and the authors should change the paper to more clearly state the change in strategy. Simply put, the new algorithm is NOT “unconstrained”.
We agree with the reviewer that our reformulation is not strictly unconstrained but bound-constrained with the bound being that parameters should be positive. We will update the title of the manuscript accordingly and will add a paragraph to clarify the terminology and better present the non-linear inequality constraints used in the original solid solution phase formalism of e.g., Holland et al. (2018).
The central problem is simply stated: minimize the Gibbs free energy of a phase as a function of the fractions of cations X_i, I=1,n on one or more crystallographic sites. The minimization is subject to two constraints:
1) The cation fractions must sum to 1:
sum_i^n X_i = 1
and
2) they must be non-negative
X_i >= 0
Constraint (1) can be implemented either by performing a constrained minimization for which (1) forms an auxiliary statement to the minimization problem OR by minimizing over the null space of the constraint (1). The latter approach is appealing because it reduces the dimensionality of the problem (by the number of crystallographic sites), and avoids the auxiliary statement of constraint.
We agree with the general description of the problem. However, the solution phase formulation as described by Holland et al. (2018) and all previous publications of the THERMOCALC group are given in terms of compositional variables, which within bounds, can results in site fraction < 0.0. There is thus a need to reformulate the solution phase using mixing site constraints. As discussed by the reviewer below, for other (simpler) thermodynamic databases, this may not be necessary.
The primary focus of this paper is on the implementation of the null space approach. And it is on the basis of the removal of the auxiliary statement (1) that they describe their modified algorithm as “unconstrained”.
However, it is NOT unconstrained. The reason is constraint (2). This must still be dealt with: the minimization must still be subject to the bounds described by constraint 2. Thus the minimization is still constrained, not unconstrained.
We agree that the problem is still bound-constrained and will clarify it throughout the manuscript.
In other words, they have removed an equality constraint, but NOT the inequality constraint.
We appreciate the opportunity to clarify our intent. Our intention was not to suggest the removal of the bound constraints of the minimization problem.
Our contribution is that we reformulate the inequality constraints of the site fraction, which are expressed as functions of compositional variables (e.g., Holland et al., 2018), into equality constraints expressed as functions of site fraction variables. Subsequently, we eliminate these equality constraints using a nullspace approach, leaving only the bound constraints on the site fraction variables.
This approach is adopted because the site fraction formulation presented in Holland et al. (2018) incorporates inequality constraints that restrict the hypercube defined by the 'compositional variables.' In other words, the hypercube defined by the bound constraints of the compositional variables does not fully represent the feasible domain. For instance, certain combinations of compositional variables can result in negative site fractions (even within bounds), thus necessitating the use of these additional inequality constraints. The focus of the paper is the elimination of these inequality constraints.
From the reviewer remarks we understand that this was insufficiently clear in the text, so we will update the manuscript accordingly in the revised version.
Of course the authors DO end up applying the inequality constraint, although this application is somewhat buried in the details (Eq. 32).
Each dimension of the problem is indeed bounded, but as described above, we no longer must invoke a solver with (internal) inequality constraints.
My argument is not with the method itself, nor primarily with the claims of greater efficiency (although see below). It is with what I think is a misleading description of the modification. I implore them to characterize their modification more clearly and correctly.
We hope that the revised version clarifies this better.
Some further comments.
1) Readers may get the mistaken impression that the null space approach is new in the context of petrological Gibbs free energy minimization codes. It is not, and an example of another minimization code that uses the null space approach is HeFESTo (Stixrude and Lithgow-Bertelloni, 2011). Nor do the authors of HeFESTo claim priority for the null space approach, but it is the closest implementation known to this reviewer in terms of intended application (i.e. to petrology).
We agree with the reviewer that the nullspace approach is not a novel concept in the context of petrological Gibbs free energy minimization and we did not intent to present it as such. Indeed, this method forms the foundation of the solution phase formalism in THERMOCALC, particularly in the formulation of “compositional variables”, which have been used since at least White et al. (2000). These variables are used to parametrize mixing site charge neutrality, enforce the condition that the sum of site fractions equals 1, and ensure that the sum of endmember proportions also equals 1.
2) The claimed superiority of their new implementation to SLSQP is puzzling to this reviewer. One reason for the source of puzzlement is that HeFESTo uses SLSQP and finds a very high (essentially perfect) rate of success and precision as documented in Stixrude and Lithgow-Bertelloni (2021).
This is likely because the problem solved using SLSQP in HeFESTo is only bound-constrained, while the one solved in MAGEMin also includes internal inequality constraints. This is why we wanted to reformulate the problem in a similar manner as in Stixrude and Lithgow-Bertelloni (2021). We will clarify this in the revised version of the manuscript.
Perhaps the reason is that HeFESTo uses the null space approach for the equality constraints and relies on the constraint facility of SLSQP only for the inequality constraint. Maybe in the current paper, the authors are instead relying on SLSQP to take care of the equality AND inequality constraints.
We think that there is a misunderstanding coming from having a different solution phase formalism in mind here.
While the solution phase formalism described in Stixrude and Lithgow-Bertelloni (e.g., 2021, as well as earlier references) is indeed constrained solely by bounds (referred to by the reviewer as inequality constraints), the site fraction formulation presented in Holland et al. (2018) is natively different. It incorporates non-linear inequality constraints that internally restrict the hypercube defined by the compositional variables bounds. For example, certain combinations of compositional variables (within bounds) can still result in a negative site fraction. Therefore, additional non-linear inequality constraints are required. Our manuscript shows how to take those into account before passing the problem to the innermost solver.
Therefore, to avoid confusion in the literature, I suggest the following test:
Perform the SLSQP minimization(s) again, but by combining SLSQP with the null space approach.
The results of the SLSQP algorithm presented in the manuscript already used the nullspace approach as the solution phase formalism presented in Holland et al., 2018 ensures that the mixing site charge neutrality, the sum (site fraction) = 1 and the sum of endmember proportion = 1 by parameterizing the solution space.
Moreover, using the SLSQP approach on bound-constrained only problems reduces to using a bounded BFGS approach.
3) A final comment on SLSQP. The authors state that SLSQP sometimes fails because it violates inequality constraints and then cannot return to the feasible space. But this should not be true. According to NLOPT documentation SLSQP is guaranteed to respect inequality constraints at all intermediate steps of the minimization.
We respectfully think that the reviewer is mixing bound-constraints and (non-)linear inequality constraints. While bound constraints are indeed always respected, this is not necessarily true for non-linear inequality constraints.
The NLopt documentation indicates the following (https://nlopt.readthedocs.io/en/latest/NLopt_Tutorial/):
“In principle, we don't need the bound constraint x2≥0, since the nonlinear constraints already imply a positive-x2 feasible region. However, NLopt doesn't guarantee that, on the way to finding the optimum, it won't violate the nonlinear constraints at some intermediate steps, while it does guarantee that all intermediate steps will satisfy the bound constraints. So, we will explicitly impose x2≥0 in order to ensure that the √x2 in our objective is real.”
So, in other words, the nonlinear constraints are not always guaranteed to be satisfied.
Perhaps the problem is that SLSQP sometimes does venture into the space where one or more X_i are exactly zero. If this is the case, it is a problem easily solved, by setting the inequality constraint instead to:
X_i >= epsilon
very much like their own implementation of the inequality constraint.
As pointed out by the reviewer, having a bound-constraint for a site fraction exactly equal to zero would pose a problem (as log(0.0) is undefined for the configurational entropy term) and we indeed added a small epsilon in the code to avoid this (see https://github.com/ComputationalThermodynamics/SandBox/blob/main/GradientBasedMinimizers/unconstrained_CG_BFGS/functions/gradient_method.jl L99 and L237-245).
We will further clarify this in the text.
Citation: https://doi.org/10.5194/gmd-2024-197-AC1
-
AC1: 'Reply on RC1', Nicolas Riel, 27 Jan 2025
reply
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 198 | 41 | 13 | 252 | 7 | 10 |
- HTML: 198
- PDF: 41
- XML: 13
- Total: 252
- BibTeX: 7
- EndNote: 10
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1