CALPHAD: Computer Coupling of Phase Diagrams and Thermochemistry 34 (2010) 452–455
Contents lists available at ScienceDirect
CALPHAD: Computer Coupling of Phase Diagrams and Thermochemistry journal homepage: www.elsevier.com/locate/calphad
On choosing a reference surface for a two-sublattice model Dmitri V. Malakhov Department of Materials Science and Engineering, McMaster University, JHE 357B, 1280 Main Street West, Hamilton, Ontario, Canada, L8S 4L7
article
info
Article history: Received 12 July 2010 Received in revised form 27 August 2010 Accepted 28 August 2010 Available online 17 September 2010 Keywords: Sublattice model Compound energy formalism Reference surface Ill-posed problem Regularization
abstract Hillert’s choice of a reference surface for the (A, B)a (C, D)c solution has vital advantages: the corresponding contribution to the Gibbs energy of a phase is straightforward to calculate; the surface is geometrically simple; it is continuously differentiable; the method can easily be generalized for multicomponent phases with any number of sublattices regardless of whether a component resides in one sublattice only or may inhabit several sublattices. Elegance, convenience, simplicity and naturalness of his construction may mask the fact that the choice of the reference frame is not unique. From the utilitarian angle, a discussion of the non-uniqueness undertaken here may be of little use, because Hillert’s method will forever be used for building the reference surfaces for phases with two or more sublattices. A sheer curiosity rather than pragmatism motivated revisiting a problem of defining the reference surface and considering Tikhonov regularization and the least-norm criterion as unconventional ways of selecting it. Mathematical particularities and peculiarities associated with these two approaches are discussed in the present work. © 2010 Elsevier Ltd. All rights reserved.
1. Introduction
all four components:
Let us consider a quaternary solution made of the components A, B, C and D, and assume that the two-sublattice model (A, B)a (C, D)c is employed for describing its thermodynamic properties. It is convenient1 to choose the sizes of the sublattices so that a + c = 1; yet, in this work a = c = 1 is assumed. The assumption does not jeopardize a generality of subsequent derivations, but makes them less cumbersome. If y and z are used for designating the site fractions of B and D, then (A1−y By )(C1−z Dz ) can be interpreted as the model of the phase written in a format conforming to the spirit of the compound energy formalism [1,2]. It can also be viewed as the phase itself having a particular composition. For a given composition, i.e., for the certain values of y and z, there are at least two ways to make the quaternary solution. First, it can be fabricated from the pure components:
(1 − y)A + yB + (1 − z )C + zD → (A1−y By )(C1−z Dz ). Second, the pure components can be used to form four binary compounds (end-members) AC, AD, BC and BD from which the phase is then obtained according to the reaction:
α AC + β AD + γ BC + δ BD → (A1−y By )(C1−z Dz ).
(1)
At first glance, the stoichiometric coefficients α , β , γ and δ in (1) can be derived from the conditions of material balance written for
1 0 1 0
1 0 0 1
0 1 1 0
0 α 1−y 1 β y = . 0 γ 1 − z 1 δ z
0364-5916/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.calphad.2010.08.003
(2)
However, the stoichiometric coefficients cannot be found, because the rank of the matrix in the LHS of (2) is equal to 3:
(1 1 0 0) − (1 0 1 0) + (0 0 1 1) = (0 1 0 1) . 1st row
3rd row
2nd row
4th row
Since (2) is an underdetermined system of equations, α , β , γ and δ cannot be uniquely defined, which means that a reference surface associated with the Gibbs energies of the end-members G0AC , G0AD , G0BC and G0BD cannot be unambiguously constructed without invoking a postulate. Hillert sagely suggested finding such stoichiometric coefficients that would provide the simplest geometrical shape of the reference surface [3]:
α = (1 − y)(1 − z ), β = (1 − y)z , γ = y(1 − z ), δ = yz .
(3)
The substitution of (3) in (2) yields an identity. The reference surface defined as Gref = (1 − y)(1 − z )G0AC + (1 − y)zG0AD
+ y(1 − z )G0BC + yzG0BD E-mail address:
[email protected]. 1 The reason behind this expediency is revealed in the Appendix.
(4)
has numerous benefits. First, it is easy to calculate; second, Gref is continuously differentiable with respect to y and z, and the
D.V. Malakhov / CALPHAD: Computer Coupling of Phase Diagrams and Thermochemistry 34 (2010) 452–455
453
derivatives are easy to compute as well; third, a generalization of (4) to the case of multicomponent phases with several sublattices is straightforward. It is clear that Hillert’s choice will never be seriously challenged, and that it will forever be used in the compound energy formalism. Obviously, there is no utilitarian reason justifying a discussion on his selection of the reference surface. Let us, however, mute the voice of pragmatism and explore alternative ways for constructing the reference surface. 2. The solution resulting from Tikhonov regularization The underdeterminacy of the system (2) is not unusual; illdefined problems are frequently encountered in the practice of scientific and engineering calculations [4]. A routine way of dealing with such problems is to use Tikhonov regularization [5]. Let us recall the essence of his method by considering an overdetermined system of linear equations Ax = y, where A is an [m × n] fundamental matrix,2 x is an [n × 1] vector of unknown parameters, and y is an [m × 1] vector of experimental observations; n is equal to the number of basis functions used. The least-squares technique is customarily employed to find a particular solution xˆ minimizing ‖Ax − y‖2 , which is the sum of squares of residuals.3 It can be shown that in this case, xˆ = (AT A)−1 AT y. Apparently, this approach does not work if A is singular. The method may also fail even if from a strict mathematical viewpoint, the fundamental matrix has a complete column rank. The reason is that the columns may be linearly dependent from the computational viewpoint, because a finite number of bytes are used for representing machine numbers. An ill-conditioned matrix results in an unstable numerical solution of the system Ax = y. In order to bypass this obstacle, it was suggested to minimize ‖Ax − y‖2 + ‖0x‖2 , where 0 is the so-called Tikhonov matrix. It can be shown that in this case, the solution is given by xˆ = (AT A + 0T 0)−1 AT y. In the simplest case, the [n × n] identity matrix I can be used as 0, but typically an effect of the Tikhonov summand is controlled by a scaling of I using a non-negative regularization parameter τ : 0 = τ I. One can make τ so large that ‖τ Ix‖2 = τ 2 ‖x‖2 would dominate and xˆ = 0 would be the solution of the optimization problem ‖x‖2 → min, to which the original problem was degraded by choosing a huge τ . This stable minimum norm solution is indeed of little interest. The case τ = 0 is of no interest as well, because, as already known, it yields an unstable solution. Fortunately, there is a range of regularization parameters within which the term ‖τ Ix‖2 ensures that a stable numerical solution exists without nullifying the presence of the term ‖Ax − y‖2 related to the actual problem in hand. Let us show how the interval of suitable regularization parameters can be found when Tikhonov regularization is utilized for solving the ill-posed problem (2). For illustrative purposes, let us take y = 0.35 and z = 0.25. For this case, Hillert’s choice (designated by the superscript H) yields α H = 0.4875, β H = 0.1625, H γ H = 0.2625 and δ = 0.0875. The Euclidean norm of this solution is 0.48752 + 0.16252 + 0.26252 + 0.08752 ≈ 0.584. If regularization is employed for solving the system of equations (2), then a particular choice of τ will result in a parameter-specific solution, which is not what one wants. Then how can we choose a suitable τ ? An inspection of Fig. 1 answers this question. If τ → ∞, then one ends up with the least-norm solution α = β = γ = δ = 0 whose substitution in (2) results in different LHS and RHS, which is hardly unexpected. If the regularization parameter decreases, then the solution first undergoes drastic changes, but then it stabilizes and asymptotically approaches α = 0.45, β = 0.2,
2 [m × n] means that the matrix has m rows and n columns. 3 ‖p‖ symbolizes pT p.
Fig. 1. Influence of the regularization parameter on the solution of the underdetermined system (2) and its Euclidean norm, and a comparison of the solution and the norm with those resulting from Hillert’s choice.
γ = 0.3 and δ = 0.05 whose substitution in (2) gives an identical LHS and RHS. In fact, there is a range of regularization parameters within which the solution depends on τ so weakly that the α(τ ), β(τ ), γ (τ ) and δ(τ ) functions are virtually constant. If τ decreases further, then sooner or later the solution will become unstable because the effect of regularization would be lost. Fig. 1 is representative in the sense that the plateau at which the solution of a system of linear equations with an ill-conditioned matrix is practically independent of τ is almost always observed when a scaled identity matrix is used. The Euclidean norm of the solution
0.452 + 0.22 + 0.32 + 0.052 ≈ 0.579
is less than that of Hillert’s solution, which is not surprising, because Tikhonov regularization favors small-norm solutions [5]. 3. Minimum norm solution If it is intended to find the least-norm solution, then there is no need in using regularization, because such a solution can be directly found by solving the following quadratic programming problem [6]:
α 2 + β 2 + γ 2 + δ 2 → min 1 1 0 0 α 1−y 0 0 1 1 β y 1 0 1 0 γ = 1 − z 0 1 0 1 δ z α ≥ 0, β ≥ 0, γ ≥ 0, δ ≥ 0.
(5)
For y = 0.35 and z = 0.25, the solution of (5) is α = 0.45, β = 0.2, γ = 0.3 and δ = 0.05. It coincides with that resulting from Tikhonov regularization, but does not require finding the plateau within which the solution is essentially invariant of τ . The least-norm solution has an interesting characteristic illustrated in Fig. 2. If the composition of the solution expressed in y and z falls within the gray right-angled triangle touching the AC endmember, then δ = 0, which means that the BD compound is not needed for forming the phase having this particular composition. A similar conclusion can be arrived at with respect to the other three gray right-angled triangles depicted, because within each of them, one of the stoichiometric coefficients is equal to zero. Intuitively, such a result is expected, because whatever the site fractions, only three end-members are needed for synthesizing the
454
D.V. Malakhov / CALPHAD: Computer Coupling of Phase Diagrams and Thermochemistry 34 (2010) 452–455
in Fig. 2 do not overlap, which means that they correspond to compositions for which a ‘‘redundant compound’’ can unambiguously be identified. Although the reference surface resulting from the least-norm solution of the underdetermined system (2) demonstrates interesting features, such a choice will never be used in practice. First, in contrast to (4), the reference surface is not given as an explicit mathematical expression (a formula), but as a numerical solution of the quadratic problem (5). Second, and more important, the surface is not differentiable. The easiest way to demonstrate this is to consider a particular case. Let us use, for the sake of determinacy, the following Gibbs energies of end-members (in Joules per formula unit): G0AC = 5000, G0AD = 12,000, G0BC = −3000 and G0BD = 7000. Then, let us construct two reference surfaces: one with the stoichiometric coefficients conforming to Hillert’s choice, and another with the least-norm coefficients. Fig. 4(a) is a graphical illustration that (4) is a continuously differentiable function. In contrast, the contour lines in Fig. 4(b) have singularities positioned along the boundaries (shown as dashed straight lines) separating the region where all stoichiometric coefficients are positive from those within which one of them is equal to zero. The existence of break points means that the reference surface is not smoothly differentiable. This, in turn, is an ultimate argument against using the least-norm approach in practice.
Fig. 2. If the composition of the system belongs to a gray right-angled triangle, then one of the stoichiometric coefficients found by means of the least-norm technique is equal to zero.
phase with a corresponding composition. Let us take a particular composition y and z corresponding to a thick dot in Fig. 3(a). Since the dot is distant from the BD corner, there is no need in using this compound:
4. Concluding remarks The location of a reference surface for a phase with multiple sublattices is not uniquely defined. In order to fix its position, an additional postulate must be invoked. Hillert’s idea to use the geometrically simplest surface was superb, because the resulting reference surface possessed numerous advantages. Although his practical choice will always be employed in thermodynamic modeling, it was interesting to revisit the problem of nonuniqueness and explore the properties of a reference surface constructed on the basis of a different postulate. Since the problem is ill-posed, it was tackled by using Tikhonov regularization. As expected, there was a range of the regularization parameter within which the solution was almost independent of it. Exactly the same solution was obtained when it was postulated that the solution had to have the minimum Euclidean norm. It was found that a direct application of the least-norm criterion led to particular
(1 − y − z )AC + zAD + yBC → (A1−y By )(C1−z Dz ). The gray right-angled triangle AC − AD − BC in Fig. 3(a) is the geometric place of points for which the conditions 1 − y − z ≥ 0, z ≥ 0 and y ≥ 0 are satisfied. Since the dot is distant from the AD corner as well, the use of this end-member is avoidable, too:
(1 − y)AC + (y − z )BC + zBD → (A1−y By )(C1−z Dz ). The inequalities 1 − y ≥ 0, y − z ≥ 0 and z ≥ 0 are fulfilled within the gray right-angled triangle AC − BD − BC in Fig. 3(b). Since the thick dot belongs to both triangles, an attempt to find the stoichiometric coefficients by minimizing the minimum one does not result in a unique solution. In contrast to this, the gray triangles
a
b
Fig. 3. Ambiguity associated with the optimization problem min(α, β, γ , δ) → min.
D.V. Malakhov / CALPHAD: Computer Coupling of Phase Diagrams and Thermochemistry 34 (2010) 452–455
a
455
z
z
b
y
y
Fig. 4. A continuously differentiable reference surface resulting from Hillert’s choice (a), and a non-smooth reference surface corresponding to the least-norm approach (b).
combinations of the site fractions for which one end-member was not needed. Ostensibly, that could be considered as an advantage, because an evaluation of the Gibbs energy of a corresponding compound could be avoided in the course of thermodynamic assessment. However, the least-norm criterion resulted in a reference surface that was not continuously differentiable, which negated an apparent benefit of having a reduced number of endmembers. Although it is clearly understood that this short communication is not very useful from the utilitarian perspective, it is hoped that the discussion of mathematical aspects related to the nonuniqueness of the choice of the reference surface would be of some interest to the CALPHAD community.
One concludes that the total number of moles of all components is equal to the total number of moles of all phases. Keeping this in mind, let us consider the following reaction:
Acknowledgement
summation of the conditions of mass balance nA = Na(1 − y′ ), nB = Nay′ , nC = Nc (1 − y′′ ) and nD = Ncy′′ yields nA + nB + nC + nD = N (a + c ). The total number of moles of the components is equal to the number of moles of (A1−y′ By′ )a (C1−y′′ Dy′′ )c only if a + c = 1.
Financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) is gratefully acknowledged. Appendix Let us consider a multicomponent system, which was formed by taking nA moles of the component A, nB moles of the component B and so on. Let us assume that these components were distributed φ among the phases α, β, . . . , that the system is composed of. If xi is the molar fraction of the component i in the phase φ , then the following conditions must be valid: ni =
−
φ
nφ x i ,
i = A, B, . . .
φ
where nφ is the number of moles of the phase φ . Let us calculate the total number of moles of all components in the system:
− i
ni =
−− i
φ
φ
nφ x i =
− φ
nφ
− i
φ
xi =
− φ
nφ .
(1mole of A) + (1mole of B) → N moles of the compound, in which xA = xB = 1/2.
(6)
It is clear that N = nA + nB = 2 and that, therefore, the compound should be denoted as A1/2 B1/2 rather than AB. Subsequently, the reaction (6) should be written as A + B → 2A1/2 B1/2 rather than A + B → AB. In the case of the reaction nA A + nB B + nC C + nD D → N (A1−y′ By′ )a (C1−y′′ Dy′′ )c
References [1] J.-O. Andersson, A.F. Guillermet, M. Hillert, B. Jansson, B. Sundman, A compound-energy model of ordering in a phase with sites of different coordination numbers, Acta Metallurgica 34 (1986) 437–445. [2] M. Hillert, The compound energy formalism, Journal of Alloys and Compounds 320 (2001) 161–176. [3] M. Hillert, L.-I. Staffansson, The regular solution model for stoichiometric phases and ionic melts, Acta Chemica Scandinavica 24 (1970) 3618–3626. [4] G.E. Forsythe, M.A. Malcolm, C.B. Moler, Computer Methods for Mathematical Computations, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1977. [5] A.N. Tikhonov, V.Y. Arsenin, Solutions of Ill-Posed Problems, Winston & Sons, Washington, 1977. [6] D. Goldfarb, A. Idnani, A numerically stable dual method for solving strictly convex quadratic programs, Mathematical Programming 27 (1983) 1–33.