Prediction of temperature-dependent yield strength of refractory high entropy alloy based on stacking integrated framework

Database establishment

This study focuses on RHEA systems, which consist of nine high melting point elements (Zr, Ti, Hf, Nb, Mo, V, Ta, W and Cr) and one light element (Al). Data on the temperature-dependent yield strength and compositions of these alloys were compiled from multiple published studies ^[7,23–26]. Data from multiple sources, preparation methods, and testing conditions were collected to enhance both the diversity and reliability of the dataset. The data were selected to encompass a wide range of chemical compositions by including alloys with varying elemental compositions and ratios. To capture the microstructural diversity of RHEAs, alloys with different grain sizes, morphologies, distributions, and phase compositions were also considered. The data were preprocessed to eliminate duplicate values, outliers, and entries that were broken within the elastic strain range. After preprocessing, the dataset consisted of 275 data points: 54 ternary alloys, 168 quaternary alloys, and 53 quinary alloys. Given the empirical criterion for testing data selection, the dataset is separated into sets for training and testing in a 4:1 ratio.

The selected feature is closely linked to the target attribute, ensuring that the model's predicted values are close to the real values within an acceptable error range. According to Ref. ^[21], the mechanical properties of high entropy alloys are highly dependent on their phases, such as hardness, yield strength, ultimate tensile strength, and so on. Hou et al. studied the relationship between the solid solution phase and valence electron concentration (VEC), mixing entropy ($$ \rm \Delta {S_{mix}} $$), mixing enthalpy($$ \rm \Delta {H_{mix}} $$), atomic radius difference(δ), electronegative difference ($$ \rm \Delta χ $$), $$ \rm \Omega $$ and other parameters of high entropy alloy, and established a model to solve the phase classification problem^[27]. The bulk modulus (K) and shear modulus (G) are directly connected to the alloy's mechanical characteristics ^[28]. K calculates the fracture resistance and demonstrates the degree to which the alloy's volume changes in response to external forces. G measures the extent of shear deformation and evaluates the propensity of fracture resistance to rise following the commencement of plastic deformation. Moreover, increasing the experimental temperature enhances atomic diffusion capacity, increases the number of vacancies in the material, and alters the grain boundary slip system, leading to significant fluctuations in strength. As a result, it is considered an important feature^[29].

Supplementary Table 1 summarizes the 21 descriptors developed in the current study. The properties of each element were taken from the appropriate sources ^[2], and the enthalpy of mixing between the two elements was determined from a table in the literature ^[30]. Because there is no evident distinction between "solvent" and "solute" in HEAs, the average weighting of elemental characteristics is used to define this random solid solution, as given in

where the alloying element proportion is expressed by $$ c_{i} $$ and the related properties by $$ P_{i} $$. Considering the tendency of some features, such as VEC, to exhibit strong localization, the "concentration mean difference" is calculated here to represent the combined effect of elemental property mismatching, as given in

where $$ \bar P $$ is the mean value.

In the process of model training, it is necessary to convert the original data into dimensionless index evaluation values. This is because large numbers can cause numerical issues, and when the model uses gradient descent, elliptical contours result in more iterations. Therefore, some algorithms require scaling of eigenvalues, such as support vector machine, integrated learning, K-nearest neighbors (KNN), etc. We use Z-score feature normalization to standardize the data:

where $$ x $$ is the original value of the data, $$ x' $$ is the normalized feature value, $$ {\rm{\mathsf{μ}}} $$ is the mean of the population data, and $$ {\rm{\mathsf{σ}}} $$ is the standard deviation of the population data.

Feature screening

Although an excessive number of distinctive features offers a physical foundation for yield strength prediction, feature duplication reduces physical insight into the problem and raises the likelihood of overfitting. As a result, feature screening is required to identify the most representative feature subset, thereby improving the model's interpretability to some extent ^[31,32].

Pearson correlation coefficient and feature importance

The presence of a large number of highly correlated characteristics causes multicollinearity, which means that they are assigned greater weights, passively decreasing the information contained in other features and so compromising the regression model's prediction accuracy. The Pearson correlation coefficient (PCC) value measures the degree of correlation between any two features and ranges from -1 to 1. The closer the PCC value is to 1, the higher the similarity between the features. The PCC value is calculated by:

(4) $$ {\rm PCC}= \frac{{\sum\limits_{i = 1}^n {({x_i} - \overline x )({y_i} - \overline y )} }}{{\sqrt {\sum\limits_{i = 1}^n {{{({x_i} - \overline x )}^2}} } \sqrt {\sum\limits_{i = 1}^n {{{({y_i} - \overline y )}^2}} } }} $$

where the numerator is the covariance of features, and the denominator is the product of the standard deviations of the two features. If the absolute value of the correlation coefficient is greater than 0.9, it implies that the two characteristics have a strong association ^[33]. Figure 1A illustrates the PCC heat map, highlighting three sets of correlated features: [("VEC", "K"), ("$$ \rm {T_m} $$", "$$ {\rm{\mathsf{ρ}}} $$", "$$ \rm {Tb} $$", "AN") and ("G", "$$ {\rm{\mathsf{σ}}} \rm G $$")].

Figure 1. Feature selection steps: (A) heat map of Pearson's correlation coefficient; (B) lollipop plot of feature importance; (C) RFE process with AIC and R as evaluation indexes. RFE: Recursive feature elimination; AIC: Akaike information content.

By calculating the PCC values, we have taken into account the relationship between the variables. The next step is to fit a predictive model using the dataset, establish the logical relationship between features and the target attribute, and measure the relative contribution of each input feature using feature importance. This process helps effectively remove variables with low values. Figure 1B visualizes the importance of each feature, with values decreasing from left to right, where perpendiculars of the same color represent combinations of features with high correlation. Features with low importance ("VEC", "$$ \rm {T_m} $$", "$$ {\rm{\mathsf{σ}}} \rm G $$", "$$ \rm {Tb} $$", "AN") were removed. Finally, 16 relatively independent features with great influence on the target were selected from the original 21 features by PCC value and feature importance.

Model construction

To select a model with strong fit and generalization performance, we split the original dataset into a training set and a testing set in a 4:1 ratio, using key features as input. We then applied several classic ML algorithms including extreme gradient boosting regressor (XGBR), extremely randomized trees (ERT), RF, adaptive boosting (AdaBoost), GPR, KNN, and LR to the training set for repeated model building. Preliminary modeling, summarized in Supplementary Table 2, indicated that the XGBR model provided superior predictive accuracy, suggesting that gradient boosting algorithms are well-suited for predicting temperature-dependent yield strength. Consequently, we included additional models in the gradient boosting framework, such as histogram-based gradient boosting regressor (HistGBR) and LightGBM regressor (LGBMR).

While traditional ML models effectively handle structured data, their noise resistance and extrapolation capabilities are often limited. The stacking EL algorithm, based on the meta-learner concept, addresses this by combining multiple ML algorithms over two or more stages to produce a model with lower classification or regression errors. This progressive, optimal learning approach has been effective in various applications, such as predicting glass formation in amorphous metals ^[36], detecting sensor faults ^[37], and supporting hydropower emergency responses^[38]. Therefore, we leverage the EL model to improve predictive accuracy in this study.

HEAs have complex compositions, and the interactions among the various elements can significantly influence the alloy's properties. One key advantage of the stacked integration framework is its ability to leverage the strengths of multiple meta-learners, enabling it to capture diverse features in the data. This approach allows for a more comprehensive consideration of factors such as alloy composition, structure, temperature, and other variables that affect yield strength, thereby enhancing prediction accuracy. In contrast, traditional methods often rely on empirical formulas derived from limited experimental data, which may not be applicable across the full range of alloy compositions and temperature conditions. Additionally, physical models may be overly simplified and fail to accurately capture the intricate interactions and microstructural changes within the alloys. Furthermore, studies ^[39–41] also demonstrated that the stacked integration framework yielded significant improvements in predicting the properties of RHEAs, offering a more reliable tool for the design and optimization of these materials.

The EL procedure involves using the highest-performing model as the base learner, training each independently, and then stacking their outputs to form a new dataset with dimensions ($$ a $$, $$ b $$), where $$ a $$ is the sample count and $$ b $$ the number of base learners. The meta-learner then fits this transformed dataset to produce the final prediction. To prevent overfitting and manage multicollinearity, we selected simple models for the meta-learner, such as Lasso or linear regression, with regularization techniques, which have proven effective across fields, including sparse robust estimation ^[42], multi-core learning prediction ^[43], and feature selection and clustering ^[44]. Figure 2 shows the flow diagram of the integrated learning algorithm.

Figure 2. The flow of the stacking algorithm.

The use of stacking EL algorithm, which integrates multiple models, can increase prediction complexity and reduce transparency. In practice, when the output deviates from expected results, it can be difficult to isolate which specific model or input variable is responsible, especially when considering the complex interactions between various alloying elements and phases in RHEAs. This lack of interpretability can undermine the model's usefulness in explaining specific material behaviors. For instance, a prediction of yield strength under extreme temperatures might be affected by phase instability or grain boundary weakening, which may not be adequately captured by the integrated model.

In ML, overfitting where the model performs exceptionally well on training data but poorly on unseen data can hinder generalization. To mitigate this, we use $$ k $$-fold CV to evaluate the model's generalization capacity, tuning parameters via grid search. Here, $$ k $$ is set to 10, as this value minimizes bias in error estimation ^[45]. Repeating 10-fold CV ten times further stabilizes our findings and improves reliability. We selected the coefficient of determination (R$$ ^{2} $$), root mean square error (RMSE), and mean relative error (MRE) as evaluation metrics. MRE, an L1-norm loss, encourages feature sparsity, while RMSE, an L2-norm loss sensitive to outliers, promotes feature selection to enhance generalization. Together, MRE and RMSE complement each other in balancing model performance.

Model comparison

Hyperparameter optimization plays a critical role in improving model performance. In this study, grid search was employed to fine-tune the hyperparameters of each ML model. For instance, with the XGBR model, hyperparameters A, B, and C were explored across a grid of 2, 640 possible combinations. Figure 3 illustrates the final XGBR model in three dimensions, displaying max depth, learning rate, and n estimators. Supplementary Table 3 lists the optimized hyperparameters for each model.

Figure 3. Grid search results for three hyperparameters of XGBR model. XGBR: Extreme gradient boosting regressor.

Model performance was evaluated using 10-fold CV on the training set, with metrics including the R$$ ^{2} $$, RMSE, and MRE. Figure 4A shows the results, where RMSE is presented in its negative form, and MRE is expressed as 1-MRE for consistency. Coverage areas represent each method, where a larger area indicates better model performance. HistGBR and XGBR exhibited superior coverage areas, encapsulating the performance of other models. ERT stood out for its strong RMSE values. Consequently, HistGBR, XGBR, and ERT were selected for further investigation.

Figure 4. Model validation: (A) radar plot for CV; (B) scatter plot of model fit on the training and testing sets; (C) test set relative error distribution diagram. CV: Cross-validation.

Following the classical ML model analysis, we examined model performance using EL. A stacking model was constructed by pairing the three selected methods and applying the lasso algorithm as a meta-learner for ten times 10-fold CV on the training data [Supplementary Figure 1]. Compared to individual ML models, EL models showed higher R$$ ^2 $$ and 1-MRE values, lower RMSE values, and shorter error bars, indicating superior accuracy and reliability. Specifically, the ERT and HistGBR combination produced a stacking model with smaller error bars for R$$ ^2 $$, 1-MRE, and RMSE, demonstrating more consistent performance and making it the final model selection. The model's predictive performance was then evaluated on the test set [Figure 4B]. Most data points are closely aligned with the diagonal, and the R$$ ^2 $$ on the test set is 0.873, demonstrating strong agreement between predicted and experimental values. R$$ ^2 $$ value of 0.87 indicates that 87% of the variance in the material properties is explained by the model, suggesting a strong correlation between predicted and observed data. Statistically, this suggests that 87% of the variance in the dependent variable (e.g., material properties such as yield strength, hardness, etc.) can be explained by the model's independent variables.

Higher R$$ ^2 $$ values, in practice, can help mitigate issues such as material wastage, production delays, and safety incidents caused by inaccurate predictions. They can also be leveraged to optimize designs and processes, enhancing performance and reducing costs. A high-precision predictive model enables engineers to more accurately assess and optimize performance, effectively screen and refine the composition and heat treatment processes of RHEAs, and ultimately advance the use of these alloys in high-temperature and more demanding environments.

Model evaluation

Model interpretability analysis

ML models, often viewed as black boxes, can obscure the mechanisms behind their predictions, limiting insights into the mechanical characteristics of RHEAs. To address this, our study uses interpretability strategies to analyze the effects of various attributes on prediction outcomes, categorizing methods as global or local based on their focus. Global methods, including partial dependence plots (PDP) ^[51] and ALE ^[52], show the average influence of features across the data distribution. Local methods, such as SHAP ^[53–55] and individual conditional expectation (ICE)^[52], focus on individual instances.

Here, we used the pdpbox package to generate PDP and ICE plots, with PDP showing the average dependency between yield strength and selected features, while ICE plots illustrate individual sample effects, minimizing confounding factors. Supplementary Figure 4 shows that all of the ICE curves are consolidated at a single starting point, making it easier to understand the connection between local and global effects. By analyzing the features from both perspectives, it becomes clear that temperature has a negative impact on yield strength, and the influence of each feature on individual samples generally follows the overall average trend. As temperature increases, the material's lattice structure changes, and atomic thermal vibrations intensify, leading to a negative correlation between temperature and the target properties. This, in turn, enhances the material's plastic deformation ability while reducing its yield strength.

To further investigate, SHAP values and ALE plots are employed, with the average effect line in the ALE plots expected to lie within the gray confidence interval at a 95% confidence level. Figure 6A shows the SHAP summary plots, which rank the importance of the model's seven attributes. The temperature distribution across the samples is highly spread, indicating that this feature has the most significant influence on the prediction. The interaction between the lattice constant a and other features has no notable effect, which is attributed to its low importance. Therefore, taking temperature, atomic size difference δ and bulk modulus K as research parameters, Figure 6B illustrates that both from a single point perspective and an average effect standpoint, the atomic size difference δ and bulk modulus K exhibit a critical point, marking a clear boundary between two distinct groups in the plots.

Figure 5. Fitted scatterplot for experimental data by the model based on ERT and HistGBR. ERT: Extremely randomized trees; HistGBR: histogram-based gradient boosting regressor.

Figure 6. Model validation: (A) SHAP value distributions for key features of the stacking model. The data points represent individual alloy samples, with each point color-coded according to the magnitude of the relevant feature. Positive SHAP values indicate that the feature increases the yield strength, while negative SHAP values suggest a decrease in yield strength due to the feature; (B) SHAP dependence scatter plot and ALE plot for single-feature δ and K; (C) 3D partial dependence graph of interaction feature [K, G]. SHAP: SHapley Additive exPlanations; ALE: accumulated local effect.

By analyzing the intersection of the interval defined by this boundary, we find that $$ δ>0.049 $$ positively influences yield strength, while $$ δ<0.049 $$ has a negative effect. The conclusion aligns with the findings of Wang et al., which shows when δ increases, HEAs are more likely to form a high-strength BCC/B2 phase^[36]. BCC-type HEAs generally have a high yield strength. Similarly, the SHAP plot for the volume modulus K indicates that when K $$ > $$ 150 GPa, it contributes positively to yield strength, and vice versa, which supports the findings of Lee et al. that an increase in the alloy's K value is often associated with improved fracture resistance during plastic deformation^[56]. Furthermore, when a univariate effect is not significant, we explore the impact of bivariate interactions on the target property using a PDP three-dimensional diagram. As shown in Figure 6C, high values of both bulk modulus K and shear modulus G have a beneficial effect on yield strength. This finding is consistent with previous studies ^[57] that a combination of high bulk and shear modulus increases the yield strength for MoNbTaTi refractory complex concentrated alloys. The combination of a high bulk modulus and high shear modulus can enhance the yield strength of RHEAs through several mechanisms, including improved cohesion, increased resistance to dislocation slip, enhanced deformation resistance, and greater stability of the alloy's crystal structure.

Yield strength optimization

Optimizing RHEA compositions for high yield strength at elevated temperatures is critical for material development. This study combined the stacking model with the dung beetle optimizer (DBO), a recent swarm intelligence algorithm by Xue et al.^[58]. DBO thoroughly explores the solution space to avoid local optima, showing robust performance. As confirmed by Giles et al., yield strength is temperature-dependent, making independent searches for optimal compositions at 1, 000 and 1, 200 ℃ effective^[59].

Table 1 presents optimized alloy compositions. For CrNbMoV and AlCrNbMoV systems at 1, 000 ℃, which exhibit baseline strengths of 1, 036 and 1, 085 MPa, respectively, optimization yielded Al$$ _{0.050} $$Cr$$ _{0.330} $$Nb$$ _{0.060} $$Mo$$ _{0.337} $$V$$ _{0.223} $$ with a strength of 1, 300.638 MPa, representing a 264.638 MPa increase, primarily due to a rise in Mo content. Similarly, the optimized composition Cr$$ _{0.275} $$Nb$$ _{0.215} $$Mo$$ _{0.349} $$V$$ _{0.161} $$ achieved 1, 242.259 MPa, a 206.259 MPa again. At 1, 200 ℃, optimized alloys within Cr-Ta-Ti-W-V, Cr-Ta-W-V systems yielded Cr$$ _{0.181} $$Ta$$ _{0.190} $$Ti$$ _{0.050} $$W$$ _{0.301} $$V$$ _{0.278} $$ and Cr$$ _{0.256} $$Ta$$ _{0.238} $$W$$ _{0.24} $$V$$ _{0.266} $$, showing yield strength improvements of 18.16% and 3.84%, compared to the original data.

Effect of element content on yield strength

After multiple DBO cycles scanning the composition space, a dataset was generated, detailing the compositions and corresponding yield strengths of various RHEAs. This dataset facilitates analysis of the relationship between element ratios and performance. Given the complexity of RHEA systems, each alloy system was analyzed individually, with results for different temperatures and systems combined into four datasets. To reduce model error impact, we sorted the data by yield strength, selecting the top 200 entries exceeding the original alloy's target. Alloy compositions were rounded to two decimal places to avoid duplication.

Figure 7A illustrates the relationship between elements and yield strength for the Al-Cr-Nb-Mo-V alloy at 1, 000 ℃. Two distinct categories emerge: Mo concentration ranges from 17% to 35% for high yield strength, with Al(at%) = 5%, Cr(at%) = 26%-33%, and Nb(at%) = 5%-9%. V shows a broader distribution without a specific range. Thus, high Cr and low Al and Nb concentrations favor stronger alloys in the Al-Cr-Nb-Mo-V system at 1, 000 ℃. Similarly, Figure 7B shows the CrTaTiWV system at 1, 200 ℃, where high performance is achieved with Cr(at%) = 15%-23%, Ta(at%) = 16%-21%, Ti(at%) = 5%, W(at%) = 26%-35%, and V(at%) = 17%-35%. Supplementary Figure 5 provides parallel coordinates for the CrNbMoV system at 1, 000 ℃ and the CrTaWV system at 1, 200 ℃, revealing similar ranges that favor yield strength.

Figure 7. Parallel coordinates plots show the values of the elements of (A) AlCrNbMoV and (B) CrTaTiWV system; Fitted scatter plot of (C) linear regression; (D) quadratic polynomial regression.

Additionally, we examined the overall impact of nine RHEA elements. Using a dataset of 276 experimental samples, linear and quadratic polynomial regressions were applied, with temperature included as a feature. The dataset was split 4:1 into training and testing sets. Figure 7C displays the linear regression model, which achieved an R$$ ^2 $$ score of 0.665 on the test set, indicating a clear relationship between yield strength and elemental content and temperature. This linear model offers a basic understanding of each element's effect. In contrast, the quadratic polynomial model achieved a better R$$ ^2 $$ score of 0.807 [Figure 7D], highlighting its superior predictive ability and suitability for quantitative yield strength predictions. Coefficients and intercepts for the quadratic terms are detailed in Supplementary Table 4.

Authors' contributions

Investigation, visualization, writing - review and editing: Yu L

Data curation, visualization, writing - original draft: Zhai J

Visualization, writing editing: Cao W

Methodology, visualization, supervision, writing - review and editing, funding acquisition: Ren J

Availability of data and materials

The detailed materials and methods used in the experiment are provided in the Supplementary Materials. Other raw data that support the findings are available from the corresponding author upon reasonable request.

Financial support and sponsorship

The authors would like to acknowledge the financial support from the NSFC (Nos. U23A2065, 52071298), the Natural Science Foundation of Henan Province (No. 232300420346), Key Research Programs of Higher Education Institutions in Henan Province (25A110001), and the Training Plan for Young Backbone Teachers of Henan University of Technology.

Conflicts of interest

Ren J is the guest editor of the Special Issue, while the other authors have declared that they have no conflicts of interest.

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Copyright

© The Author(s) 2024.