An Interpretable Boosting-based Predictive Model for Transformation Temperatures of Shape Memory Alloys

1. Introduction

Shape memory alloys (SMAs) exhibit reversible martensitic transformations from a low-temperature, low-symmetry martensite phase to a high-temperature, high-symmetry austenite phase [1]. As a result of this reversible transition, they exhibit several behaviors, including the Shape Memory Effect (SME) [2,3], superelasticity [1,3], elasto-/magneto-caloric effects [4,5,6,7,8,9,10,11], tunable thermal expansion [12,13], stiffness [14,15], magnetic exchange bias [16], and spin glass transitions [17,18,19].

Advances in SMAs are intermittent. A significant discovery has only occurred once a decade since the late 1960s. For example, the discovery of NiTi in the 60s [20,21], Cu-based SMAs in the 1970s [22], superelasticity in the 1980s [23,24], magnetic SMAs in the 1990s [25], $\beta$-Ti SMAs in the 2000s [26,27,28,29,30,31], and multi-component high-temperature SMAs in the 2010s [29]. There is a need to develop a methodology that will accelerate the discovery of new SMAs, thereby reducing the discovery cycle from decades to just a few months.

Conventional alloy discovery has relied on a tradeoff system where the gain in one target leads to another worsening. Due to the numerous tradeoffs that often must be employed in such systems, optimizing multiple characteristics becomes even more challenging. It primarily targets binary systems where interactions between elements are often known or at least predictable. In conventional alloys, the quest for finding tradeoffs to tune specific characteristics has motivated the exploration of the multi-principal element, or high entropy alloy design space [32,33,34,35,36]. A similar paradigm has emerged in functional alloys with several compositionally complex SMAs, and Multi-Principal Element Multi-Functional Alloys (MPEMFA), which may demonstrate superior performance compared to less complex counterparts [37]. As an example of these quests, recently discovered quaternary and quinary MPEMFAs exhibit reversible thermally induced martensite transformation, SME, and superelasticity above 700°C, the highest temperature superelasticity reported to date [38].

Efficient discovery of new compositionally complex alloys is challenging due to the vastness of the chemical and microstructural space and extreme sensitivity to minute changes in chemistry and/or microstructure [39]. The successful discovery of MPEMFAs is at the expense of exploring a very high-dimensional, non-linear, and highly singular mathematical space, in which querying a single point experimentally comes at a high cost.

Expanding the chemical discovery space may lead to new discoveries, but it may also result in undesirable properties. This is especially true when considering their solid-state actuation properties and martensitic transformation temperature characteristics for optimizing actuation properties under transformation temperature constraints. As an example, in TiPt and Pd systems, higher alloy complexity led to smaller irrecoverable/permanent strains and actuation thermal cycling, but at the expense of a larger thermal hysteresis [40]. As chemical complexity increases, atomic interactions make the understanding of the new design space more complicated as numerous factors may be able to yield the same result.

Current alloy development frameworks are insufficient, as simulation-assisted frameworks can be computationally expensive, do not readily incorporate data from experiments, and are sequential. Traditional high throughput (HTP) combinatorial approaches are incapable of dealing with the high dimensionality and complexity of the MPEFMA space. In addition, HTP approaches are ‘one-shot’ or ‘open loop’ schemes without an iterative framework, and are suboptimal in resource allocation. State-of-the-art machine learning (ML)/AI-driven approaches may be unsuitable in cases where the materials space is highly non-linear, sparsely queried, or very singular. The use of modeling, combined with high throughput experiments, has the potential to vastly increase the rate of discovery in these new systems due to the newfound ability to apply active learning methods to find trends in small elemental variations, which could never be observed before.

Avoiding computationally expensive calculations such as first-principles calculations [41,42], or theoretical studies [43,44,45], multiple research groups have tried to use classical machine learning (ML) models and, more recently deep neural networks to create surrogate models, connecting chemistry to transformation temperatures. Most of the models are fitted to small training data sets that are less diverse in terms of the included elements and different SMA systems, limiting the application space of these models [45,46,47,48,49,50,51]. On the other hand, some data sets are enriched with computational data stemming from simulation models such as CALPHAD [52], which limit the results to the uncertainties and assumptions of the underlying calculations [53,54].

Although the importance of processing parameters on martensitic transformation is known, only limited studies [46,47,48,55] have considered the effect of processing parameters in their model to avoid the complexities added by these features. Wang et. al. in their recent work [56] attempted to cover the shortcomings of the literature and have combined composition-based featurization with a fractional encoding of the materials. However, the complexity of their model makes tracking and providing feature importance calculations a cumbersome task. In addition, their model limits the feature space to chemistry-related information, and processing parameters are not included. As a result, having a model that covers several alloy systems and is based on a relatively diverse data set, which considers the features containing the composition-based information and the manufacturing parameters, is missing in the literature.

Considering the aforementioned missing pieces to have a complete and comprehensive model, we are structuring our work to build a surrogate model that can fit the multi-dimensionality of a diverse yet limited data set. The precise feature engineering enabled the possibility of making predictions for various alloy systems. In addition, the generated features and Shapley values provided the chance to shed light on the physics of the transformation.

Our work is structured as follows; We first describe the data set, the background of the ML algorithm, and the surrogate ML model in the methods section. Subsequently, we will discuss our results and the novelty of this work providing physics-related information about the features affecting the transformation temperatures.

Figure 1. The number of data points available for each alloy system in the database. Please note that each alloy system is composed of multiple compositions, and each composition might have multiple data points with different processing parameters. The figure was plotted using Matplotlib [57].

Figure 1. The number of data points available for each alloy system in the database. Please note that each alloy system is composed of multiple compositions, and each composition might have multiple data points with different processing parameters. The figure was plotted using Matplotlib [57].

2. Methods

In this section, we first elaborate on the structure of the data set under the study, data pre-processing, and feature generation procedures to make our reported results reproducible. More importantly, we introduce the basics of the adopted ML surrogate model based on CatBoost, for which we will demonstrate its generalizable prediction performance with respect to the corresponding evaluation metrics in our experiments in Section 3. Last but not least, by adopting the Shapley value [58], we also provide the feature importance analysis to interpret the learned CatBoost model for SMAs transformation temperatures.

2.1. Data Preparation

A database with information on composition, processing parameters, and thermodynamic responses was constructed by combining the data mined from the literature, and our in-house experimental results [59]. Composition data was represented in atomic percentages at the resolution of 0.1at%. Processing parameters were used, including the parameters in heat treatment, rolling, and extrusion. Both composition and processing parameters construct an initial feature list for the ML surrogate model. Thermodynamic responses were limited to transformation temperatures as these are the targets for this study.

After cleaning the data (i.e., removing data points with duplicated and NaN values) and removing inconsistent values (i.e., removing data points where atomic percentages exceeded 100, Martensite Start (Ms) > Austenite Finish (Af), Martensite Finish (Mf) > Austenite Start (As), As > Af, and Mf > Ms), 2,494 data points for 237 different shape memory alloy were acquired. These data points have 13 different processing parameters. The dataset in the compositional space has data concentration around specific compositions. For some compositions, we had as low as one experimental measurement data point per composition; for others, there were more than ten data points per composition. This fact creates an under-sampling of the composition space for some regions. Still, it doesn’t imply an imbalance in the dataset, as the data points with the same composition have different processing parameters. To add another layer of model evaluation for generalizability, we divided the data set into the train/test set and the validations set. The train/test set contains compositions with at least two experiments per composition and is divided into training and test sets for appropriate training of the model and to avoid overfitting the training data. On the other hand, the validation set is the set of compositions with only one experiment. The set is used to evaluate the model further. The number of samples available for alloy systems in each data set is shown in Figure 1.

We employed three featurization methods for the present work: processing-related, composition-related, and alloy-system-related. Compositions were featurized with the Jarvis database [60] using the CBFV [56,61,62,63,64,65] Python [66] library. We added these features based on previous experiments and density functional theory, ML, and classical force-field calculations, as a result, increasing the total number of features to 3,079. For more detailed explanations of the features generated in this paper, please refer to JARVIS-Tools [67] and CBFV [68] documentation.

We note that our final feature set also contains eight categorical features to help our ML surrogate model better predict thermodynamic properties for different SMA systems. These categorical features are NiTi-based, NiTiPd-based, NiTiCu-based, NiTiHf-based, NiTiNb-based, NiTiZr-based, NiTiB-based, and alloy systems. This set of categorical features provides information on SMA systems with their main element combinations. As an example, the alloy system for Cu2Ni3Pd5Ti10 is NiTiPdCu. Other categorical features are subcategories of alloy systems. In the case of Cu2Ni3Pd5Ti10, it belongs to the NiTi-based, NiTiPd-based, and NiTiCu-based categories, so the value for these features is True, and the value for the four remaining categorical features (i.e. NiTiHf, NiTiNb, NiTiZr, and NiTiB-based) is False. As a result of feature engineering, the database contained 3,087 features ready for training.

2.2. CatBoost Surrogate Model

The number of alloy systems, significant change in thermodynamic targets with slight composition variation, number of features, and presence of categorical and numerical features in the database indicate a complex non-linear learning space. Decision-tree-based ML models are good candidates for capturing nonlinear composition\processing\structure property-space mapping relationships in a high-dimensional learning space. Different boosting and bagging strategies, built upon tree-based weak learners, such as Random Forests [69], XGBoost [70], LightGBM [71], and CatBoost [72], have been developed to achieve state-of-the-art prediction performance. Convenient in working with categorical data and having multiple targets are the other reasons that we chose CatBoost as our final ML surrogate model. In the following subsections, we introduce the background of the ML algorithm, model evaluation methods, and hyperparameter tuning by Bayesian optimization (BO).

2.2.1. Background of the ML Algorithm

Boosting [73] is a family of ML algorithms that combines multiple weak learners to build a strong learner for prediction. Decision trees [74], which recursively partition the input space and learn the decision rules for each subspace based on training data, are usually chosen as weak learners or base learners because of their lightweight nature and low computational costs.

Gradient boosting [75], as a member of the boosting family, utilizes the gradient of a pre-defined loss function to combine the weak learners. Unlike other implementations of gradient boosting, including XGBoost [70] and LightGBM [71], CatBoost [72] tries to fight potential prediction shifts with an ordered boosting method and has been reported to achieve more reliable prediction performance.

In CatBoost, the strong learner $M_t$ is updated as [72]:

$$M_{t+1}\left(\mathbf{x}\right)=M_t\left(\mathbf{x}\right)-\alpha h_t^{*}\left(\mathbf{x}\right),$$

where $h_t^{*}\left(\mathbf{x}\right)\in {arg\; min}_{h\left(\mathbf{x}\right)}\left(E_{\mathbf{x}}\left[{(h\left(\mathbf{x}\right)-\frac{\partial L(f,\mathbf{x},y)}{\partial f}{|}_{f=M_t})}^2\right]\right)$ is the estimation of $\frac{\partial L(f,\mathbf{x},y)}{\partial f}{|}_{f=M_t}$ with the weak learner and $\frac{\partial L(f,\mathbf{x},y)}{\partial f}{|}_{f=M_t}$ is the gradient of the loss function w.r.t the current strong learner $M_t$,

$\alpha$ is the learning rate to control the update step size. $L(f,\mathbf{x},y)=\mathcal{L}(f,\mathbf{x},y)+J\left(f\right)$, which is the loss function connecting the prediction to the target, usually contains an L2 regularization term $J\left(f\right)$ to prevent from overfitting.

As another way to resist prediction shift and speed up the training procedure, Bayesian bootstrap [76] is applied in CatBoost in each iteration on the training dataset to sub-sample and re-weight the data and balance model inputs.

2.2.2. Loss Function and Evaluation

The choice of the loss function $\mathcal{L}(f,\mathbf{x},y)$ is critical in the training process to achieve accurate predictions. To measure the prediction accuracy for a single target, Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and coefficient of determination (R2 score) are broadly used. The R2 score, measuring the dependency relationship between the prediction and target, is defined as [77]:

$$R^2=1-\frac{{\sum }_i{(y_i-M\left({\mathbf{x}}_i\right))}^2}{{\sum }_i{(y_i-\overline{y})}^2},$$

where $\overline{y}$ is the average over all the target values.

RMSE measures the difference between the target and the prediction by [78]:

$$\mathrm{RMSE}(M\left(\mathbf{x}\right),y)=\sqrt{\frac{{\sum }_i{(M\left({\mathbf{x}}_i\right)-y_i)}^2}{N}},$$

where $M\left({\mathbf{x}}_i\right)$ is the prediction of the i-th sample, $y_i$ is the corresponding target, and N is the number of samples. MAE, as another broadly used metric, is defined as [78]:

$$\mathrm{MAE}(M\left(\mathbf{x}\right),y)=\frac{{\sum }_i|M\left({\mathbf{x}}_i\right)-y_i|}{N}.$$

A perfect predictor (always correct) would have a 0 RMSE and MAE with an R2 score of 1.

When it comes to multiple target prediction, we can adopt RMSE to MultiRMSE as follows:

$$\mathrm{MultiRMSE}(M\left(\mathbf{x}\right),y)=\sqrt{\frac{{\sum }_i{\sum }_j{(M{\left({\mathbf{x}}_i\right)}_j-y_{i,j})}^2}{N}},$$

where $M{\left({\mathbf{x}}_i\right)}_j$ is the prediction of the j-th target of the i-th sample (${\mathbf{x}}_i$), $y_{i,j}$ is the corresponding target, and N is the number of samples.

2.2.3. Hyperparameter Tuning

Fine-tuning the model’s hyperparameters can be a pathway to pursuing higher prediction accuracy. Before training the model, we tuned the hyperparameters of the CatBoost algorithm using the BO approach [79], utilizing the Bayesian Optimization [80] Python library. The optimization process was performed on the training set, and the negative square root of the mean_squared_error was set as the target of the BO algorithm to optimize. The resulting parameters used to build the ML model are illustrated in Table 1.

By considering the hyperparameters as the input to a black-box function whose output is the corresponding model performance, the task of fine-tuning the hyperparameters can be viewed as a black-box optimization problem. BO [79] that iteratively queries new input and eventually optimizes the prediction performance is suitable under this scenario. BO models the black-box performance function as a Gaussian process (GP) $\mathcal{GP}\left(\mu \right(·),Cov(·,·\left)\right)$. An acquisition function guides the selection of new query points based on the updated posterior of the GP condition on the current observations.

The choice of acquisition function is critical and affects the optimization efficiency. Expected Improvement (EI), Entropy Search (ES), and Upper Confidence Bound (UCB) are examples of well-studied acquisition functions [79]. We focus on the UCB acquisition because of its low computational cost and strong performance, which takes the form of:

$$\alpha \left(x\right)=\mu \left(x\right|\mathcal{D})+\beta \sqrt{Cov(x,x|\mathcal{D})},$$

where $\mathcal{D}$ denotes the set of current observations, $\mu \left(x\right|\mathcal{D})$ and $Cov(x,x|\mathcal{D})$ represent the mean and variance of the updated posterior GP, respectively. UCB balances the exploration and exploitation by the hyperparameter $\beta$, where higher $\beta$ values lead to a higher level of exploration while lower $\beta$ lead to exploitation more. In our experiments, $\beta =2.576$.

2.3. SHAP-Based Model Interpretation

Besides prediction accuracy, model interpretability is another goal when developing AI/ML surrogate models, particularly for scientific machine learning (SciML), as in this paper. To better understand the composition-processing-property relationships, we also calculate the SHapley Additive exPlanation (SHAP) [58] for each feature, enabling information-theory-based feature prioritization. The corresponding SHAP values capture the importance of the related features in determining materials’ thermodynamic responses, which provides insight into SMA physics.

The SHAP Value $\varphi (k,M)$ for the feature k in a model M is defined as:

$$\varphi (k,M)=\sum _{S\subset F/k}{w}_{S,F}[M^{S\cup \left\{k\right\}}\left({\mathbf{x}}^{S\cup \left\{k\right\}}\right)-M^S\left({\mathbf{x}}^S\right)],$$

where F is the set of all the features, $M^S\left({\mathbf{x}}^S\right)$ is the retrained model with only the features in a subset S instead of all the features.

The retraining of $2^{\left|F\right|-1}$ models to calculate the SHAP value of each feature, however, is computationally infeasible in most applications. The original paper [58] proposed several approximation methods to tackle this exponential computational complexity challenge with significantly reduced retraining time. SHAP Value can be viewed as a weighted average of the difference between the retrained model with and without feature k for each subset S, providing a quantitative feature importance metric capturing its contribution to the final model prediction.

3. Results and Discussion

As detailed in Section 2.1, we have collected 2,494 data points for a set of 237 SMA compositions composed of 26 different alloy systems. The training was performed on 80% of the training/test data set, which contains 137 SMAs from NiMnCoIn, NiTi, NiTiAu, NiTiCu, NiTiCuAu, NiTiCuCo, NiTiCuNb, NiTiCuZr, NiTiHf, NiTiHfCu, NiTiHfFeCu, NiTiHfZr, NiTiHfZrCu, NiTiNb, NiTiPd, NiTiPdB, NiTiPdCo, NiTiPdCu, NiTiPdHf, NiTiPdSc, NiTiPdTa, NiTiPt, NiTiZr, NiTiZrAl, and NiTiZrB alloy systems, and the performance of the model was evaluated on the test set (20% of the above data set). Furthermore, a “validation” set that contains 100 SMAs from NiTi, NiTiCu, NiTiCuZr, NiTiHf, NiTiHfZr, NiTiNb, NiTiPd, NiTiPdCu, NiTiPdSc, NiTiPdTa, NiTiPt, NiTiZr, and NiTiNbZr alloy systems was held out, excluding from the above training/test data set to better evaluate the generalizability of the trained CatBoost surrogate.

Figure 2. Predicted and actual values of Ms, Mf, As, and Af temperatures for the test set. The figure was plotted using Matplotlib and seaborn [81]. The visually distinct colors generated by iWantHue [82].

Figure 2. Predicted and actual values of Ms, Mf, As, and Af temperatures for the test set. The figure was plotted using Matplotlib and seaborn [81]. The visually distinct colors generated by iWantHue [82].

The average $R^2$ score, MAE, and RMSE of the model predictions were adopted for performance evaluation, calculated using scikit-learn [83]r2_score, mean_absolute_error, and mean_squared_error functions. On the training set (80% of the training/test data set), the average $R^2$, MAE, and RMSE values of Ms, Mf, Af, and As temperature predictions are 0.95, 18.47°C, and 29.97°C. The same evaluation metrics for the test set (20% of the training/test data set) are 0.95, 18.75°C, and 26.78°C. The high $R^2$ score for the training and test set

is an indication of not having under/overfitting for the trained model. Figure 2 shows the model’s performance for the test set. The $R^2$ scores for Ms, Mf, Af, and As temperatures are 0.95, 0.96, 0.96, and 0.95 respectively. The MAE for the same targets is 17.39°C, 18.32°C, 18.00°C, and 21.29°C. The RMSE for the same transformation temperatures are 25.10°C, 25.23°C, 25.97°C, and 30.45°C. Also, when we look into more detailed thermodynamic predictions for different alloy systems, it is clear that our CatBoost model with customized input features not only correctly predicted transformation temperatures of NiTi-based SMAs but also did a good job predicting NiMn-based alloys.

We have further tested the model’s performance on compositions unknown to the model to investigate how well this model can predict the transformation temperature of alloys, which the model has never seen during the previously described training/test procedure. This analysis is performed on the data set we referred to as the “validation” set with the distribution shown in Figure 1. This data set contains 100 NiTi-based shape memory alloys, and the performance of the model on this data set is shown in Figure 3. These data points are unique from the train/test data points either in the sense that they have a unique composition or processing parameters.

The comparison plot of the transformation temperatures for this validation set also shows the proximity of the predicted and true transformation temperatures. The $R^2$ score, MAE, and RMSE of the model for Ms temperature are 0.92, 18.42°C, and 22.96°C, for Mf temperature are 0.89, 22.84°C, and 30.53°C, for Af temperature are 0.92, 21.33°C, and 27.44°C, and for As temperature are 0.89, 26.76°C, and 34.22°C.

It is worth mentioning that although the model was never trained with any NiTiNbZr alloy system, it could predict the targets for these alloy systems. Also, only 6 out of 26 alloy systems had more than 150 data points in the training/test set, where the number of data points for the NiTiCu alloy system is lower than that of the validation set. Another interesting fact, as shown in Figure 4, is that the model could detect the sequence in transformation temperatures without defining any constraints. Accordingly, the model could understand and keep the following sequence: Ms < Af, Mf < As, As < Af, and Mf < Ms.

Last but not least, the interpretation of the SHAP-based model was performed. The absolute Shapley values for each feature are averaged to create the SHAP feature importance plot, shown in Figure 5, which depicts the magnitude of feature attributions. This plot suggests that for the NiTiX alloy system, transformation temperatures are mainly a function of chemical and electronic properties such as atomic mass (atom_mass), atomic radii (atom_rad), molar volume (mol_vol), heat of fusion (hfus), melting point (mp), boiling point (bp), electronegativity (X), electron affinity (elec_aff), and polarizability (polzbl). Moreover, the final heat treatment temperature (FinalHT_Temp) as a processing parameter is among the top 10 most important features affecting this transformation temperature.

The uneven amount of SHAP values for targets (consider the width of the colorful area for each feature) suggests that not all features evenly affect transformation temperatures, which shows the complexity of the model in making predictions.

To better illustrate how changing values of a feature affect transformation temperature predictions, the summary plot for Ms temperature is depicted in Figure 6 as an example. It shows that high values of average atomic mass and electronegativity difference (avg_atom_mass_subs_x), population standard deviation of heat of fusion and electronegativity difference (dev_hfus_subs_X), molar volume and heat of fusion difference (avg_mol_vol_subs_hfus), population standard deviation of atomic radii and atomic mass quotient (dev_atom_rad_divi_atom_mass) increase the $M_s$ temperature. On the contrary, low values of population standard deviation of electronegativity and electron affinity difference (dev_X_subs_elec_aff) and summation of polarizability and atomic radii (sum_polzbl_add_atom_rad) tend to increase this target. In the case of processing parameters, the high values of final heat treatment temperature can increase or decrease the martensitic temperature, based on its contrary effects on the microstructure and precipitates of each composition [84]. Hence, it offers a handy tool for selecting alloying elements for designing SMAs with high transformation temperatures.

When it comes to the computational costs, training the model for 2,394 data points with 3,087 features and 4 targets on a Windows laptop equipped with Intel^® Core^TM i7-11800H processor and 16 GB Random Access Memory (RAM) took 3 minutes and 21 seconds. The model proves to be fast, accurate, and simple - it has competitive performance and can be an excellent choice for time-sensitive studies.

4. Conclusions

Employing advanced ML techniques and domain knowledge, we were able to build a model to predict transformation temperatures of SMAs, provide insight into the important features contributing to each prediction, and shed light on the physics of thermodynamic transformation of SMAs. The results of the present work are as follows:

• Proper feature engineering with all processing parameters, composition-based features, and alloy system categorical features is the key to building a robust model that is able to predict transformation temperatures of a vast range of SMAs with a relatively low amount of data. In our study, the model could predict the transformation temperatures of 100 SMAs in the validation set after training by 2,394 data points for 137 alloys in the train/test set.
• The results shows that CatBoost is a promising ML surrogate model for the thermodynamic prediction of complex alloy systems. Considering its tree-based nature, it can adequately handle diverse data sets with a potentially high number of features, which is the case in most materials science studies. It also performs well with a relatively limited number of training data points. This tree-based boosting algorithm also provides a natural path to model interpretability.
• The dependency of transformation temperatures on processing parameters, chemical, and electronic properties was shown and discussed.