Researchers have used the (calculation of phase diagram) CALPHAD method to solve the forward phase stability problem of mapping from specific thermodynamic conditions (material composition, temperature, pressure, etc.) to the associated phase constitution. Recently, optimization has been used to solve the inverse problem: mapping specific phase constitutions to the thermodynamic conditions that give rise to them. These pointwise results, however, are of limited value since they do not provide information about the forces driving the point to equilibrium. In this paper, we investigate the problem of mapping a desirable region in the phase constitution space to corresponding regions in the space of thermodynamic conditions. We term this problem the generalized inverse phase stability problem (GIPSP) and model the problem as a continuous constraint satisfaction problem (CCSP). In this paper, we propose a new CCSP algorithm tailored for the GIPSP. We investigate the performance of the algorithm on Fe–Ti binary alloy system using ThermoCalc with the TCFE7 database against a related algorithm. The algorithm is able to generate solutions for this problem with high performance.
The designers of novel materials require an understanding of phase stability in order to assess the feasibility of a material and how it changes during processing. The calculation of phase diagram (CALPHAD) method has enabled researchers to develop databases that contain pertinent thermodynamic information on specific alloys and associated phases . In this method, the thermodynamics of phases are described through mathematical models fitted to experimental data.
Many calphad software packages developed to date are well-suited for solving the forward phase stability problem: calculating the phase stability state of a multicomponent, multiphase system given a set of thermodynamic conditions (processing conditions such as composition, temperature, and pressure). However, because design is inherently an inverse process, we are interested instead in the inverse problem: determining the thermodynamic conditions that give rise to a desired phase stability state.
where the set is the feasible region in the decision space (thermodynamic conditions), and f is a function (e.g., a CALPHAD model) from to that maps thermodynamic conditions to the thermodynamic degree-of-freedom of interest. The optimization problem in Eq. (1) can be solved using conventional techniques. In Ref. , mesh-adaptive direct search (MADS)  algorithms are used to find extrema in phase stability. More recently, the Arroyave group has used genetic algorithms (GAs) to solve similar problems .
Although the location of extremal points, or any isolated point(s) in the phase constitution space, is sometimes of interest to material scientists [10,11], more generalized knowledge is desirable for the materials discovery and design process. Specifically, we are interested in identifying the region or set in the design space that gives rise to a desirable range of thermodynamic conditions, we refer to this problem as the generalized inverse phase stability problem (GIPSP).
For example, a designer may want to know what thermodynamic conditions (a set of designs) result in the presence of σ phase in steels, since any amount of σ phase is undesirable due to its embrittling effects (as little as 3 wt.% reduces impact toughness by more than half ) and its drastic deterioration of the stability against corrosion. A compact representation of this set would be useful to material designers as a constraint on the search space that is easy to evaluate. In cases when uncertainty in the appropriate CALPHAD model is high, materials designers may be interested in finding a set of potentially desirable solutions that they can refine and prune through experimentation. In general, materials designers are interested in identifying the set of all the thermodynamic conditions that could produce the set of arbitrary phase constitutions. In the language of engineering design, materials designers would like to apply set-based approach  to the materials discovery process.
In the context of set-based design, Shahan and Seepersad suggested the use of Bayesian network classifiers to model the satisfactory region . In Ref. , the technique is extended to multilevel design problems (multilevel models are those that are hierarchically nested) and applied to materials design. An extension is developed by Rosen  that enables the incorporation of process knowledge. However, the focus of these works is on the design methodology concerning multiscale design rather than an efficient algorithm approach for approximating the satisfactory set, which is the aim of this paper.
We model the GIPSP as a continuous constraint satisfaction problem (CCSP) , a type of constraint satisfaction problem (CSP) where all the variable domains are continuous real intervals. The GIPSP is to map specific regions in multidimensional phase constitution spaces to ranges of values of thermodynamic conditions space. The solution set to this problem is ranges of values (enclosures) rather than discrete points. The GIPSP is particularly challenging because the search space is highly nonlinear and discontinuous, the CALPHAD model is a black-box since its details are not available, and the problems of interest may be multidimensional (>10).
Typical algorithmic approaches for CCSPs (Fig. 1), such as those based on interval arithmetic , require accessible analytical problem formulations. However, the GIPSP involves a nonanalytical CALPHAD model. Methods such as design of experiments (DOE)  or inductive design exploration (IDEM) [20,21] are not tractable for high-dimensional problems since they sample the search space. Another approach is to use adaptive sampling schemes to replace the expensive to evaluate constraint with a surrogate model [22–25]. However, these approaches have difficulty representing discontinuities in the search space and suffer from the curse of dimensionality . To address this limitation, Basudhar and Missoum developed the explicit design space decomposition (EDSD) method. The EDSD method relies on a support vector machine (SVM) classifier to model the design space constraint boundary, rather than approximate the constraint function. The EDSD method has been shown to accurately model the constraint boundary on several test problems within only a fewer iterations.
However as we will show that the SVM technique used in EDSD has difficulty representing the solution to a typical GIPSP. In a GIPSP, it is undesirable to densely sample the unsatisfactory region since it comprises the vast majority of the space. In the case where the samples are imbalanced (many more from one class than another), the SVM will be unreliable in the undersampled region. In cases with imbalanced sampling, the support vector domain description (SVDD) technique has been shown to be more appropriate . In this paper, we present an algorithm that instead relies on the support vector domain description (SVDD)  classifier, which leads to improvements in both computation time and solution quality of the GIPSP.
Generalized Inverse Phase Stability Problem as a Continuous Constraint Satisfaction Problem
where is a n-tuple of variables defined in the thermodynamic conditions space, and is the Cartesian product of the corresponding domains, where Ii is a real interval for . The set of constraints that must be satisfied is denoted as , which is a t-tuple of constraints. Each constraint Cj is a pair Rj, Sj, where Rj is a relation on the variables in . The relation Rj defines the consistent value combinations. Specifically, in the context of the GIPSP, Rj is an inequality or equality in terms of the function mapping between thermodynamic conditions and phase constitutions (i.e., constraints defined on the CALPHAD model). The set is an unordered k-tuple of distinct variables, where k is the arity of Rj. In other words, Sj is simply the tuple of variables that participate in the constraint Rj.
where the superscripts lb and ub denote the variable lower and upper bounds, respectively.
Since mapping from the thermodynamic conditions space to the phase constitution space is highly complex, discontinuous, and nonanalytical, satisfying the user-defined constraints C is nontrivial. In the motivating problem, we are interested in the set of all the solutions, , to the CCSP where all the variable domains are continuous real intervals and all the numeric relations are equalities and/or inequalities. The constraints are, in general, highly nonlinear.
Existing CSP Algorithms.
Classical methods for solving CSPs such as backtracking , iterative broadening search , and limited discrepancy search  are intended for problems with discrete domains and are not efficient for solving CCSPs. To apply these to a CCSP, one must discretize the variable space to an enumerable set . Such an approach would be intractable for most GIPSPs because of the high-dimensional (often n > 10)  of many materials design scenarios. Most techniques developed specifically for solving CCSPs are based on interval arithmetic, branch and bound, or the root inclusion test. In one of the first examples of set-based design, Ward proposed the use of interval arithmetic to search for the feasible set of designs efficiently . Subsequently, this work was extended to more general set-based representations . Hu et al. proposed a method that uses generalized interval to solve for the feasible set [35,36]. The NUMERICA  modeling language in particular guarantees correctness, completeness, and certainty. Devanathan and Ramani presented a polytope-based method, where the constraints are transformed into ternary constraints . Then, the so-called consistency method is used to prune the search space. These techniques require an analytical or closed-form expression to determine whether a subsection of the search space contains feasible solution [28,39,40]. Since the CALPHAD model is a “black-box” in the sense that the details are not accessible, methods that rely on interval arithmetic or constraint decomposition cannot be used for the GIPSP.
Design of Experiments.
Identifying the feasible set (solving a CSP) is a key challenge in many set-based design methods. Design of experiments (DOE) can be used to test the constraint model at a set of finite points in the design domain. A limitation is that the points will not lie on the boundary of the constraints. The inductive design exploration method (IDEM) addresses this limitation [20,21]. The IDEM consists of discrete point evaluations (discrete sampling) of the design process, performing inductive, top-down feasible space exploration based on metamodels. The aim is to obtain a robust solution with respect to model uncertainty. To approximate the constraint boundary, the design space is sampled using a typical DOE. The true constraint boundary will lie between the satisfactory and unsatisfactory discrete points. A root-finding technique is used to find the location of the boundary between these points. The constraint is then represented as a set of boundary points and points inside the feasible space. Although this method generates points on the boundary of the constraint, it does little to reduce the computational burden of the initial DOE, which suffers from the curse of dimensionality.
where is the joint probability density function of all the random parameters. This integral may be approximated using Monte Carlo simulation (MCS), but this approach can be computationally expensive when a large number of samples are required or the constraints are expensive to evaluate. A principle focus in RBDO is the efficient approximation of the solutions to the constraint g(x). An approach common in the literature is to replace the expensive to evaluate constraint g(x), with an inexpensive surrogate model generated from sampled data. The most straightforward methods sample the space globally, which becomes prohibitively expensive at high dimensionality. An approach for reducing the number of samples required to construct the surrogate model is adaptive sampling. Adaptive sampling approaches incrementally improve the surrogate model by sampling points that most effectively improve the model [22–24]. These approaches attempt to efficiently “train” a surrogate model over the entire search space. However, to approximate Eq. (2), it is only necessary to know where the constraint is satisfied. This insight led to the efficient global reliability analysis (EGRA) method. In EGRA, samples are generated by maximizing the expected feasibility function, which provides an indication of how well the true value of the response is expected to satisfy the constraint . Several other criteria can be found in the literature [41–43].
As argued by Basudhar and Missoum, the principle limitations surrogate-based approaches are response discontinuities and the curse of dimensionality . The discontinuities in the search space are problematic for gradient-based techniques. To address this limitation, they proposed a method referred to as explicit design space decomposition (EDSD), which does not approximate the response of the limit state function, instead constructs an explicit constraint boundary around the design variables . Thus, the EDSD method is unaffected by discontinuities and other irregularities in the search space. The EDSD method relies on the use of support vector machines to construct an explicit boundary of the constraint in the design space. Support vector machines (SVMs) are a machine-learning technique for supervised learning associated with two-class classification. That is, given a set of training examples, each example is marked as belonging to one of the two classes. The SVM uses the training examples to build a model that assigns new (unobserved) examples to one of the categories. The SVM is initially trained on a sample of points. Then, an adaptive strategy is used to generate points that likely lie on the constraint boundary.
The EDSD method has been shown to approximate the limit state function successfully for a number of test problems . However, because the EDSD method relies on the SVM method, its performance deteriorates when an SVM is not an appropriate representation of the boundary. The SVM treats evaluated samples as a two-class data set: feasible and infeasible observations. For a good performance, the SVM technique requires both classes to be well sampled [46,47].
This limitation is significant in the GIPSP, where, intuitively, few process parameter combinations will yield a desirable material. Thus, the solution to a GIPSP is typically small relative to the search space. In this scenario, it is undesirable to evenly sample both satisfactory and unsatisfactory classes, especially for high-dimensional problems where a balances sampling of the entire search space may be computationally expensive. If to reduce computational expense, we undersample the unsatisfactory space, the resulting SVM will be unreliable in the undersampled unsatisfactory space and will tend to have high rates of false positives. In cases with imbalanced sampling, the support vector domain description (SVDD) technique has been shown to be more appropriate . Whereas the SVM scheme is a two-class classifier, the SVDD scheme is a one-class classifier intended to address the case where the training data are mainly from a single category.
Support Vector Domain Description
For a detailed description of the method for formulating the Wolfe dual problem see Ref. . For each data point, xi for i = 1, 2,…, n, there are three possible classifications:
It is inside the hypersphere, which is indicated by βi = 0.
It is on the boundary of the hypersphere, which is indicated by 0 < βi < c.
It is an outlier outside of the hypersphere, which is indicated by βi = c.
A new test point, z, is inside the domain description if the distance from the feature space image of test point to the hypersphere centroid is less than the radius of the hypersphere. The expression for classification, Eq. (9), is a simple algebraic expression that is fast to evaluate. In fact for the Gaussian kernel function, the first term is equal to 1, and the last term can be precomputed since it is independent of z.
where (the index i again enumerates both target and outlier-data). See Ref.  for a detailed exposition of the SVDD method with negative examples.
To prevent weighting variables with large magnitudes more than those with lower ones in this comparison, the training data are centralizing (scale all data to a −1 to 1 range), which improves the SVDD model. An important benefit of the SVDD method is that it can be constructed incrementally and decrementally . This allows for a relatively inexpensive update procedure to be used when new members are added or removed from the SVDD. Figure 2 is an illustration of the SVDD method on a two-dimensional data set in the thermodynamic conditions space.
To underscore the difference in performance between the SVM and SVDD classifiers, we develop two test problems where one class (the satisfactory region) is (a) large and (b) small relative to the domain, see Fig. 3. In each case, we created a training data set with 50 satisfactory and 50 unsatisfactory random examples. In example (a) where the regions are proportional, the samples are relatively balanced, while in (b), they are imbalanced (the unsatisfactory region is undersampled).
We use both examples to train an SVM and SVDD classification model. We used the SVM algorithm in MATLAB's Statistics and Machine Learning Toolbox , and SVDD model in DD_Tools  developed in the Pattern Recognition Laboratory at Delft University of Technology, Delft, Netherlands. The results are illustrated in Fig. 3. In scenario (a), the SVM model outperforms the SVDD model. This is expected since the SVDD tends to generate conservative model; recall that the SVDD finds the minimum radius hypersphere around the target-data. In contrast, in scenario (b) where the data are not balanced, the SVM classifier results in a significant overestimation. In the latter case, the more conservative SVDD scheme has better performance. This basic insight motivates an adaptive sampling scheme based on SVDD rather than SVM for the GIPSP.
Our aim is to develop an adaptive sampling scheme based on the SVDD technique for approximating the constraint boundary. The basic idea of the proposed algorithm is that the true boundary of the satisfactory region is approximately parallel to our current best guess. In the proposed algorithm, we search along the direction perpendicular to the SVDD boundary (our current best guess) for a point on the true boundary of the satisfactory region using a root-finding method. Figure 4 is an illustration of this idea. An initial point is selected along the SVDD boundary along with an initial step size. The initial step size is selected using information from the SVDD. If the end point is outside of the satisfactory region, a root-finding method is used to find a point on the boundary. The proposed method is described in detail in Algorithm 1. We assume that the designer has available small number of designs that satisfy the specified phase state properties, i.e., the constraints. The initial samples may be obtained from prior equilibrium experimental data or found using conventional optimization techniques. In the case studies presented in this paper, we use random sampling to generate n data points in the design variable space. The randomly generated samples are assigned labels , such that y = 1 or y = −1 according to whether or not they satisfy the constraints, respectively. The set of indices Ib is initialized to the empty set.
where t is a dummy variable, are the training data, and for any , the set of support vectors. To prevent the case where the algorithm attempts to “reuse” a previous initial point, Ib should also include the indices corresponding to previous initial sample points. Because the Gaussian Kernel is an inexpensive function, Eq. (12) is fast to evaluate. It is important to note that the initial point xa is not guaranteed to be satisfactory. This can occur when the SVDD modeling parameter q (see Sec. 4) is set too “loose.” As a result, it is necessary to evaluate the initial point xa against the constraints to determine its label, denoted as ya.
It is desirable to choose a step size γ large enough to cross the boundary of the satisfactory region. If the label ya = 1, we should step outside of to find additional satisfactory points, “growing” the model. On the other hand, ya = −1 indicates that is optimistic, and we should step inside of to find additional unsatisfactory points, “tightening” the model. Choosing an appropriate initial step size has a significant impact on algorithmic performance. One consideration is that the initial step should not take us too far from the current SVDD, . To address this, we limit the step size to , which limits the step size according to the size and shape of along the direction of d. Another consideration is that during search, we should expect disjointed “clusters” in . In this situation, we should take a step size that is at the feature space midpoint of the clusters. If disjointed clusters exist along the direction d, their feature space midpoint is at . We use these concepts to determine the initial step size according to Algorithm 2.
The next step is to search for a boundary point using line search (spec. bisection search), see Algorithm 3 for a description. If the initial step xb did not cross the boundary of the true solution, that is, ya = yb, the bisection search algorithm terminates and returns the training set X and Y containing only the samples xa and xb. Else, the algorithm uses bisection search to reduce the size of the interval between xa and xb until it is less than some user-defined error tolerance, ε.
Finally, the training set X and Y are updated to contain the samples points that were evaluated against the constraints, i.e., points for which a label y was generated. The set of indices Ib is also updated to contain indices corresponding to the true boundary points found (within tolerance ε). This process is repeated user-defined N times, a termination rather than a convergence criteria. It would be possible to develop a convergence criteria based on the change in at each generation; this is left as future work. Without a convergence criteria, analysis of computation complexity is less meaningful but still worth considering. In the case of the GIPSP, evaluating a design against constraints (ClassLabel in Algorithms 1–3) is the elementary operation, all others are lower order. Let ε0 and ε be the maximum interval length and error tolerance for the bisection search, respectively. The time complexity is , where N is the number of iterations to be performed.
Algorithm 1 SVDD-Based Sampling Algorithm
1. procedure SVDDsample()
5. for to Ndo
6. train SVDD model with X, Y ▷ Eq. (8)
7. ▷ Eq. (12)
10. ▷ Concatenate
Algorithm 2 Take Initial Step
1. procedure TakeInitStep()
2. gradient of ▷ Eq. (13)
3. ▷ r2 from Eq. (9)
6. ifya = −1 then ▷ If xa is not satisfactory
10. return xb, yb
Algorithm 3 Bisection Search
1. procedure BisectionSearch()
2. ▷ Concatenate
3. ▷ Counter
6. ▷ Midpoint
9. ifya = ycthen ▷ Update interval
13. ▷ Update counter
14. returnX, Y, I
We evaluate the performance of the proposed algorithm on Fe–Ti binary alloy system, see phase diagram in Fig. 5. Given the thermodynamic conditions, the phase compositions are computed using ThermoCalc with the TCFE7 database. Recall that the GIPSP is the triple . The thermodynamic conditions (design variables) are , where x1 is the mass percent of Ti, x2 is the mass percent of Fe, and x3 is the temperature (Kelvin). The search domain is defined as
Results and Analysis.
The proposed algorithm is motivated by the GIPSP, which is typically multidimensional and features response discontinuities that are problematic for gradient-based search techniques. To evaluate the performance of the proposed algorithm on this class of problem, we compare to EDSD, which is intended for constraint satisfaction problems with similar characteristics. The principal difference is that the proposed algorithm is intended to address the case where the satisfactory region is undersampled (as we expect is the case with most GIPSPs).
In both algorithms, performance is dependent on the initial training data. To account for this, each case was tested for 30 trials with random initial training data. For each trial (in either test case), we randomly sampled the domain to find ten feasible and ten infeasible sites. The same data are used to initialize each algorithm. These initializing functions are not counted toward the overall function count of either algorithm.
Since both EDSD and the GIPSP algorithm use binary classifiers, we evaluate solution quality using the precision and recall metrics commonly used in pattern recognition . We generate a set X of 106 random samples in the design space domain defined by Eq. (15). We then find , the subset that satisfies the user-specified conditions, , in Eq. (16) for test case 1 and Eq. (17) for test case 2. Next, we find the subset that is classified as belonging to the satisfactory set, according to the classifier (SVM for EDSD and SVDD for the GIPSP algorithm). We compute the true positives, true negatives, false positives, and false negatives as
Lower values of misclassification rate are preferred, however, one must be cautions when interpreting this measure. For example, in a case where the satisfactory region is very small, classifying the entire region as “unsatisfactory” will have a low misclassification rate.
We calculate each performance metric at several intervals for each algorithm. We report the mean values and 95% confidence interval of 30 trials for both test cases in Figs. 6 and 7, respectively. In test case 1, the GIPSP algorithm produces a higher precision approximation for any given number of function evaluations. This is not unexpected since the SVDD is more conservative than the SVM technique. Both algorithms converge to a similar measure of recall and a low level of misclassification error. In test case 2, the GIPSP algorithm generates high-precision solutions and converges to a solution with high recall and low misclassification rate. However, the precision of the EDSD solution does not improve significantly after 100 iterations, and the recall measure becomes worse. Further, the solution quality is highly variable across each iteration, resulting in large confidence intervals.
For illustrative purposes, we include Figs. 8 and 9, which depict the progression of each algorithm at the (a) 100, (b) 250, and (c) 400, function evaluations for each test case. The shaded region represents the true solution to each CCSP found through an exhaustive search. An attempt was made to illustrate trials that are representative of the mean values in 6 and 7. However, we should not that in both test cases (but especially for test case 2), the performance of the EDSD algorithm varied significantly. Therefore, no illustration of a single result can be truly illustrative of the typical results. Taking these limitations into account, the illustrations still provide some valuable insight into the performance of each algorithm.
As can be seen in Fig. 8, the GIPSP algorithm maintains high precision during its progression, including few false positives. The EDSD, on the other hand, produces more optimistic estimations of the satisfactory region, initially. The higher precision of the GIPSP solution is reflected in the “tighter” classifier boundary. The precision at (c) 400 function evaluations for the EDSD algorithm is considerably higher than the med.
Note that Fig. 9 is focused on the solution space (which is quite small), the search space for this test problem is defined by Eq. (15) and depicted in Fig. 5. Test case 2 is intended to highlight the limitations of EDSD method on problems where the satisfactory region is small relative to the search space. In such cases, the SVM technique used in EDSD tends to overestimate the satisfactory region. In Fig. 9, the overestimation occurs near the true solution (shaded portions), however, this is not always the case. Figure 10 is an illustration of the results from another trial of the EDSD algorithm in test case 2.
Discussion and Summary
In this paper, we have presented a novel algorithm for approximating all the solutions to a CCSP with nonisolated solution where the satisfactory region is small relative to the search space. The algorithm uses the SVDD technique combined with a sampling strategy to gradually develop the solution. The motivation for the algorithm is the general inverse phase stability problem of mapping user-specified regions in multidimensional phase constitution space to ranges in values of thermodynamic conditions, which we term the GIPSP. In the GIPSP, one class (the satisfactory region) is small relative to the other (unsatisfactory region). For scalability, it is desirable to undersample the unsatisfactory region, since it comprises the vast majority of the space. This motivates the use of the SVDD method in the algorithm since it is able to more accurately (in terms of precision and recall) model scenarios with imbalanced training data.
We investigated the performance of the algorithm on Fe–Ti binary alloy system using ThermoCalc with the TCFE7 database. Using this system, we formulated two test cases. In the first test case, the solution set is nonconvex; in the second, the solution set is small relative to the search space. We compare the performance of the GIPSP algorithm to the EDSD algorithm, which uses a related classification scheme, namely, SVM. The performance of each algorithm on the test problems was measured as the precision, recall, and misclassification rate. In both test problems, the GIPSP algorithm is able to converge to a solution with high precision and recall. The EDSD algorithm, however, had significant difficulty in approximating the solution to test case 2. This is likely the result of the limitations of the SVM technique used in EDSD. The SVM technique is known to underperform in cases with imbalanced training data sets. In test case 2, the satisfactory region is small relative to the search space, resulting in an imbalanced training data set.
Future work should also investigate the performance of the algorithm on problems of higher dimensionality that are more representative of real-world materials design problems.
This work was supported by the National Science Foundation and the Air Force under Grant No. EFRI-1240483. The authors would like to thank Paul Mason from ThermoCalc for providing critical thermodynamic data, which made this research possible.
- a =
centroid of hypersphere
- b =
centroid of feature space hypersphere
- c =
- Cj =
- d =
direction perpendicular to SVDD boundary
- f =
thermodynamic conditions → phase constitution
random parameter joint pdf
- g =
- Ib =
set of indices corresponding to boundary points
- Ii =
- K =
- n =
number of data points
- N =
dimensionality of thermodynamic conditions space
- Nfn =
number of false negatives
- Nfp =
number of false positives
- Ntn =
number of true negatives
- Ntp =
number of true positives
- q =
Gaussian kernel parameter
- r =
- R =
- Rj =
relation on the variables involved in constraint Cj
- Sj =
scope of the constraint Cj
- SV =
set of support vectors
- t =
- X =
training data set
- xi =
- xi =
a vector in the design variable space
- Y =
set of training data labels
- yi =
training data label
- z =
- βi =
- γ =
- ε =
- ε0 =
maximum interval length
- ξi =
- Φ =
data space → feature space
n-tuple solution to
t-tuple of constraints
constraint satisfaction problem
n-tuple of variables