Original data consists of two aspects: one aspect is the national standard “Designation and composition of titanium and titanium alloys” (GB/T 3620.1-2016), which provides information on titanium alloy grades and their chemical composition. The other aspect involves extensive literature, where 66 relevant sources are carefully selected, primarily focusing on forged or rolled bars. The data comprises 60 titanium alloys, encompassing 18 elements, as shown in Supplementary Appendix SA1. To represent these alloys effectively, the proportion of each element is considered as a feature set consisting of 18 dimensions.

Meanwhile, heat treatment has a significant effect on the microstructure and properties of titanium alloys. In this study, two heat treatment systems are retained from selected literature, encompassing parameters such as initial heat treatment temperature and duration, subsequent heat treatment temperature and duration. The employed heat treatment system encompasses solution, aging, and all their possible combinations. Consequently, a total of four distinct dimensions are considered in this analysis. In cases where no heat treatment or only one type of heat treatment is applied, data points without any corresponding treatments are assigned a value of zero.

The original dataset comprised a total of 397 samples across 22 dimensions. In this study, the model predicts the ultimate strength as the performance indicator for titanium alloy, and the corresponding ultimate strength values were collected for each sample.

The essential properties of the elements in the alloy play a crucial role in determining its performance, thus highlighting their significance. In this study, feature engineering is employed to design the fundamental characteristics of titanium alloy, encompassing parameters such as melting point, density, atomic weight, atomic number, electronegativity, and atomic radius. However, considering that there are 18 elements in the dataset and each element corresponds to six essential characteristics, this leads to a total of 108 features or dimensions, resulting in a dimensional challenge for the dataset. To address this issue effectively while preserving relevant information integrity, a weighted summation method is adopted to reduce these 108 dimensions down to six dimensions. Consequently, each titanium alloy is associated with a set of essential features comprising weighted values for key attributes including melting point (Tm), density (ρ), atomic weight (u), atomic number (Z), electronegativity (X), and atomic radius (R). These weighted values are calculated using Formulas 1–5 through 6. Tm = ∑ i = 1 n T m i × w i (1) ρ = ∑ i = 1 n ρ i × w i (2) u = ∑ i = 1 n u i × w i (3) Z = ∑ i = 1 n Z i × w i (4) X = ∑ i = 1 n X i × w i (5) R = ∑ i = 1 n r i × w i (6) Where T m i represents the melting point value of the ith element, and w i represents the weight of the melting point for that specific element. The value denotes the content proportion of each element in titanium alloy, with a total of 18 elements (n = 18). Similar formulations are used for others. Figure 1 presents the distribution of physical property constants in titanium alloy. The values of constants are presented in Supplementary Appendix AS2.

Through the process of original data collection and feature design, the dataset consists of a total of 28 features, that is, 28-dimensional. The detailed information about these features is presented in Table 1.

In this study, feature selection encompasses both supervised and unsupervised analysis. Within the realm of unsupervised analysis, feature correlation analysis is employed to calculate the correlation coefficients between any given features X and Y, as depicted in Formula 7. r X , Y = Cov X , Y D X × D Y (7) Where Cov (X,Y) denotes the covariance between features X and Y, D(X) represents the variance of feature X, and D(Y) represents the variance of feature Y.

In the supervised analysis, the Minimum Redundancy Maximum Relevance (MRMR) algorithm is employed to compute the mutual information between the feature set and ultimate strength, quantifying feature correlation and ranking them accordingly (Peng et al., 2005). For comprehensive selection in the final feature choice, both supervised and unsupervised analysis outcomes are utilized with model accuracy as the target.

Features originate from diverse scales and encompass multiple dimensions in the dataset. To mitigate the influence of dimensionality, data necessitates processing. Data processing techniques are commonly categorized into normalization and standardization. In engineering applications, standardization is typically preferred. Within this dataset, the data without heat treatment is imputed with zeros. When normalization is performed [0,1], the feature with heat treatment would be skewed towards extreme values of 1 and 0, failing to accurately represent the overall distribution of samples in the data set. Henceforth, standardization is selected as a data processing approach, such as employing Formula 8, which centers the dataset around a mean value of 0 while maintaining a normal distribution with a standard deviation of 1. Z = C − μ σ (8)

The feature vector Z is represented as the standardized version of the original arbitrary feature vector C, where μ represents the mean value of feature vector C, and σ represents its standard deviation.

The dataset in this research is divided into three parts: the training set, the validation set, and the test set. The training set is utilized for model training, while the validation set serves the purpose of model evaluation, hyperparameter tuning, and mitigating overfitting risks. On the other hand, the test set is exclusively employed for model evaluation and testing.

Support vector machine regression (SVR) is a classic machine learning algorithm for regression proposed by Drucker et al., in 1997 (Drucker et al., 1997), and further developed by Smola et al., in 2004, where the theoretical framework and implementation of SVR was presented (Smola and Schölkopf, 2004). Similar to the SVM classification algorithm, the hyperplane used for classification often cannot completely divide the sample space, which can lead to overfitting even if it is possible to achieve complete division. In SVR, an error tolerance interval band ε is defined, and samples within this band are considered correct predictions without loss calculation. By partitioning the sample space D = {(x₁,y₁), (x₂,y₂),… (x_n,y_n)}, SVR aims to find an optimal hyperplane that minimizes the discrepancy between f(x) and y, as shown in Formula 9. f x = w T x + b (9) Where w represents the normal vector indicating the direction of the obtained hyperplane, and b denotes the displacement term representing the distance from the original point of the hyperplane, SVR computes the loss between f(x) and y outside the interval band. The regression problem is formulated in Eq. 10. min w , b 1 2 w 2 + C ∑ i = 1 n l ε f x i − y i (10)

The regularization constant C > 0 is utilized to balance the minimization of the normal vector with the minimization of the error, while l ε denotes the error tolerant interval incorporating ε insensitive loss function as defined in Formula 11. l ε z = 0 , if    z ≤ ε z − ε , otherwise (11)

By incorporating relaxation variables, the Lagrange multiplier method, and kernel functions, SVR can be mathematically formulated as Eq. 12. f x = ∑ i = 1 n α ^ i − α i k x i , x j + b (12) Where α i represents the Lagrange multiplier, the estimate of α i is denoted as α ^ i , k x i , x j denotes the kernel function, which is a fundamental technique in SVR, and refers to the transformation that maps the indivisible features of the original space to a higher-dimensional divisible space. Since the kernel function implicitly defines the feature space, the specific form of feature mapping remains unknown for a given sample space. Therefore, selecting different kernel functions often yields varying performance outcomes for regression models within specific sample spaces; an inappropriate selection may lead to reduced model performance.

The kernel function is considered as a crucial hyperparameter in this research, and the model is optimized by tuning this hyperparameter. Several commonly used kernel functions are selected, including linear, Gaussian, quadratic, and cubic kernels defined in Formulas 13–16. By adjusting these kernel functions, the optimal SVR regression model is identified. k x i , x j = x i T x j (13) k x i , x j = ⁡ exp   − x i − x j 2 2 σ 2 (14) k x i , x j = x i T x j 2 (15) k x i , x j = x i T x j 3 (16)

Gaussian process regression (GPR) is a Bayesian non-parametric probabilistic regression model that utilizes a kernel function. The application of Gaussian processes for data fitting was first introduced by O’Hagan in 1978 (OHagan, 1978). Rasmussen et al. (2005) provided a comprehensive theoretical framework for the Gaussian process regression model in machine learning in 2006. The problem of nonlinear Gaussian regression is formulated using Eq. 17. y = f x + ε (17) Where ε represents the additive noise term, ε ∼ (0, σ 2 ), where σ 2 can be estimated from the sample data. The variable y denotes the predicted value corrupted by noise, while f(x) follows a Gaussian distribution as described in Eq. 18. f x ∼ GP m x , k x , x ′ (18) Where m(x) represents the mean function and takes the value 0. The covariance function k(x,x’) captures the relationship between two features x and x’. In this regression model, Bayesian and maximum likelihood estimation methods are employed to solve for optimal parameters.

The covariance function, denoted as k(x,x’), plays a pivotal role in the model and significantly influences its performance. It is also referred to as the kernel function k(x,x’). In this study, the kernel function is treat as a hyperparameter and optimize the model by adjusting five different types of kernels: quadratic rational (Formula 19), square exponent (Formula 20), Matern 5/2 (Formula 21), Matern 3/2 (Formula 22), and exponent (Formula 23). k x , x ′ = σ f 2   1 + x − x ′ T x − x ′ 2 α l 2 − α (19)

The symbol α denotes the scaling coefficient, σ f represents the standard deviation of the sample, and l signifies the size of the feature length. k x , x ′ = σ f 2 ⁡ exp   − 1 2 x − x ′ T x − x ′ l 2 (20) k x , x ′ = σ f 2 1 + 5 x − x ′ T x − x ′ l + 5 x − x ′ T x − x ′ 3 l 2 exp   − 5 x − x ′ T x − x ′ l 21 k x , x ′ = σ f 2 1 + 3 x − x ′ T x − x ′ l exp   − 3 x − x ′ T x − x ′ l (22) k x , x ′ = σ f 2 ⁡ exp   − x − x ′ T x − x ′ l (23)

The artificial neural network is a network formed by connecting artificial neurons according to a specific topology, which originated from the MP model proposed by McCulloch et al., in 1943 (McCulloch and Pitts, 1943), as well as the perceptron neural network model introduced by Rosenblatt in 1958 (Rosenblatt, 1958). These models marked the beginning of the development boom for artificial neural networks. However, it was proven in 1969 that perceptrons were incapable of solving higher-order predicates, such as the XOR problem (Minsky et al., 1969). Consequently, the field of artificial neural networks experienced a decline until Hopfield’s proposal of the Hopfield neural network model in 1982 and Rumelhart et al.’s introduction of the BP algorithm in 1986 (Hopfield, 1982; Mcclelland et al., 1986; Rumelhart et al., 1986), which sparked a research upsurge. In this study, a backpropagation (BP) neural network is employed for regression (NNR), where each neuron’s threshold was set as θ = (θ₁, θ₂,… θ_m). Here, m represents the number of non-input layer neurons within the neural network structure. The current output Y_j of neuron j is expressed using Formula 24. Y j = f ∑ i = 1 k w ij X i − θ j (24) Where f is the activation function, i=(1,2, … k), j=(1,2, … ,m). The value of k corresponds to the number of neurons located above the current neuron j. The weight w_ij denotes the synaptic connection strength from the ith neuron in the preceding layer to the jth neuron, X_i represents the output value of the ith neuron in the preceding layer, and θ_j signifies the activation threshold of the current neuron.

The weights and thresholds of the BP neural network are iteratively calculated and updated using the gradient descent strategy, leading to the attainment of optimal model values. In this study, the neural network regression model is empirically defined as depicted in Figure 2A, comprising one input layer, one output layer, and three hidden layers. The input layer consists of 28 neurons, while the output layer comprises a single neuron; each hidden layer encompasses 10 neurons.

The neural network regression model encompasses numerous crucial hyperparameters, which are devised based on empirical knowledge that may not be deemed as ideal or relatively optimal. Hence, the optimal hyperparameters are fine-tuned by adjusting the parameters of the neural network model and iterating through Bayesian optimization to minimize mean squared error. The range for parameter adjustment is presented in Table 2. The regularization intensity value of 378 corresponds to the aggregate number of training and validation samples.

The CART algorithm, proposed by Breiman et al., in 1984 (Breiman et al., 1984), is a classic machine learning algorithm. It comprises two processes: the generation of decision trees and the pruning of decision trees. In the sample feature space, the data is divided into M units (R₁, R₂, … R_M), where each unit i has an output value defined as c_i. The regression tree model is represented by Formula 25. f x = ∑ i = 1 M c i I x ∈ R i (25) Where I(x∈R_i) represents the adaptation function for spatial feature partitioning, which is optimized by minimizing the square error and determining the optimal output value of each partition unit as the average of all observed values. The optimal feature j and segmentation point s on the partitioned feature are solved using Formula 26. min j , s min c 1 ∑ x i ∈ R 1 j , s y i − c 1 2 + min c 2 ∑ x 2 ∈ R 2 j , s y i − c 2 2 (26)

The optimal value of c₁ in R₁ is indicated by Formula 27, while the optimal value of c₂ in R₂ is denoted by Formula 28. c ^ 1 = ave y i | x i ∈ R 1 j , s (27) c ^ 2 = ave y i | x i ∈ R 2 j , s (28)

By evaluating the loss function, pruning is performed iteratively from the leaf nodes to the root node of the generated decision tree. The pruned subtree, denoted as {T₀,T₁, … T_k}, is determined based on a calculated subtree loss function using Formula 29. C α T = C T + α T (29)

The cost of pruning, denoted as C(T), is the square error associated with any subtree T, while |T| represents the complexity of the model in terms of the number of leaf nodes. Here, α is a weight parameter that balances model fitting and complexity, ultimately determining the generalization ability of the model.

In this study, the termination condition of the CART algorithm was defined as the minimum leaf size, which represents the number of samples in a leaf node. Initially, this value is set to 12 based on empirical knowledge for model regression. However, it is important to note that this initial setting may not be optimal or ideal. The hyperparameter for the minimum leaf size of the CART regression model was adjusted accordingly. Bayesian optimization and iterative processes were employed to optimize the hyperparameter based on minimizing mean squared error. Consequently, the minimum leaf size was set to [1,378/2].

The Boosting tree algorithm, proposed by Friedman et al., in 2000 (Friedman et al., 2000), is considered one of the most high-performing ensemble learning algorithms. Boosting tree is an ensemble algorithm based on decision trees and utilizes a boosting technique that employs additive models and forward distribution algorithms, following a sequential approach. The weak learners are combined into strong learners by linearly combining basis functions with weights. In this boosting tree regression model, which uses the CART algorithm as the basis function, which is defined by Formula 30. f M x = ∑ m = 1 M T x , Θ m (30) Where T(x,Θ) is the decision tree generated by CART algorithm. M represents the number of decision trees and Θ_m represents the weight of the mth decision tree. The current model of the boosting tree is defined as f_(m-1)(x), and thus, Formula 31 illustrates the weight of the mth decision tree. Θ ^ m = ⁡ arg min Θ m ∑ i = 1 N L y i , η f m − 1 x i + T x i ; Θ m (31)

The loss function is defined as the mean squared error, and the optimal weight corresponds to the weight that minimizes the mean squared error. Here, η denotes the learning rate. The boosting tree model involves several crucial hyperparameters. Based on empirical knowledge, the minimum leaf node is set to be 8, with a total of 30 learners and a learning rate of 0.1 for boosting tree regression.

Random forest, proposed by Breiman in 2001 (Breiman, 2001), is a bagging ensemble algorithm based on decision trees. Bagging involves generating multiple decision trees by randomly selecting and replacing feature and sample sets, and the predicted values of all decision trees are averaged during prediction, thereby parallelly combining the basis functions. In this random forest regression model, the CART algorithm is employed as the basis function and consider several important hyperparameters. Based on empirical knowledge, the minimum leaf node hyperparameter is set to 8 and use 30 learners for random forest regression.

However, the empirical values in the regression models of the boosting tree and random forest may not be inherently optimal or ideal. The hyperparameters of the model were adjusted and iterated using Bayesian optimization to optimize the optimal hyperparameters, based on the minimum mean squared error. The range of hyperparameter adjustments is presented in Table 2.

The root mean squared error is a widely used metric for evaluating regression models, which quantifies the average degree of deviation between predicted and true values as expressed in Formula 32. RMSE = ∑ i = 1 n y i − y ^ i 2 n (32) Where n denotes the number of samples, y_i represents the ith true value to be predicted, and y ^ i denotes the ith predicted value. The formula reveals that a smaller RMSE value indicates a closer proximity between the predicted and true values, thereby signifying an enhanced fitting degree and performance of the regression model.

The mean absolute error is a commonly employed metric for evaluating regression models, quantifying the mean discrepancy between predicted and true values as defined in Eq. 33. MAE = ∑ i = 1 n y i − y ^ i n (33)

The formula reveals that a smaller MAE value corresponds to a narrower discrepancy between the predicted and true values, indicating an enhanced level of fitting and performance for the regression model.

The coefficient of determination, a crucial evaluation index of regression models, quantifies the correlation between dependent and independent variables while representing the explanatory power of the regression model for the dependent variable. The calculation is presented in Formula 34. R 2 = 1 − SSE SST = 1 − ∑ i = 1 n y i − y ^ i 2 ∑ i = 1 n y i − μ 2 (34)

The sum of squares of error (SSE) represents the discrepancy between the predicted values and the true values in the model. On the other hand, the total sum of squares (SST) represents the overall deviation between all predicted true values and their mean value. Here, μ denotes the average value of all true values. It is evident from this formula that R ² typically ranges between 0 and 1, indicating a prediction error lower than that of mean reference. As R ² approaches 1, it signifies a higher goodness-of-fit for the regression model, implying superior performance. Conversely, when R ² is negative, it suggests poor predictive ability and greater prediction errors compared to mean reference.

The mean squared error is a widely used metric for evaluating regression models, quantifying the discrepancy between predicted and true values as computed in Formula 35. MSE = ∑ i = 1 n y i − y ^ i 2 n (35)

The formula reveals that a smaller MSE value corresponds to a closer proximity between the predicted and true values, indicating an enhanced level of fitting and performance for the regression model.

The training time reflects the model’s complexity and encompasses the time taken from the initiation to completion of training. On a consistent hardware platform, shorter training times indicate higher efficiency in model learning. For models with substantial time requirements, the training time serves as a crucial evaluation metric.

The unsupervised feature correlation coefficient results are presented in Figure 3A. It is evident that features Z and u, as well as X and Cr, exhibit a relatively high degree of correlation with correlation coefficients exceeding 0.9, which is a correlation threshold of experience.

The ranking results of feature correlation based on the supervised MRMR algorithm are presented in Figures 3A, B comprehensive analysis is conducted by considering the features associated with the supervised comparison. Specifically, Z exhibits a feature correlation score of 0.1395, u has a score of 0.1270, X is assigned a score of 0.1083, and Cr demonstrates a score of 0.1899. These features are all evaluated as significant contributors to the model’s performance. In order to construct a high-precision prediction model, all 28 features are retained without any exclusion.

After standardization, the data for each of the 28 features conforms to a normal distribution with a mean value of 0 and a standard deviation of 1, a random selection of 4 features was made to construct the histogram, as illustrated in Figure 3C.

In the support vector machine regression model, the hyperparameters were optimized to include linear, Gaussian, quadratic and cubic kernel functions. The resulting performance metrics including RMSE, MAE, R ² and MSE as well as training time T are presented in Table 3.

The optimal hyperparameter of the model is determined to be the Gaussian kernel function, as evidenced by the smallest values of RMSE, MAE, and MSE, along with an R ² value of 0.86 (Table 3). The model performance is unsatisfactory when the hyperparameter selects the cubic kernel function, as indicated by an R ² value of −148.39, which exceeds the mean reference error and indicates poor prediction accuracy. Figure 4A illustrates the comparison between real and predicted values for this particular response.

In the Gaussian process regression model, the hyperparameters of the model were adjusted to incorporate rational quadratic, squared exponential, Matern 5/2, Matern 3/2, and exponential kernel functions. The performance metrics including RMSE, MAE, R ², MSE, and training time T were presented in Table 3.

The optimal hyperparameter of the model is determined to be the exponential kernel function, as evidenced by the smallest values of RMSE, MAE, and MSE, along with an impressive R ² value of 0.95 (Table 3). The corresponding responses between real and predicted values are visually depicted in Figure 4B.

The results of the neural network regression, including RMSE, MAE, R ², MSE, and training time T, are presented in Table 4. Additionally, Figure 5A illustrates the iterative outcomes of the hyperparameter optimization algorithm. It is evident from the figure that convergence occurs during the third iteration and yields optimized hyperparameters as follows: two hidden layers with a Tanh activation function; 257 neurons in hidden layer 1 and 216 neurons in hidden layer 2; Intensity of regularization set at 0.2791. The optimized neural network regression model is depicted in Figure 2B.

The optimized NNR model has exhibited significant performance enhancement, as evident from the results presented in Table 4. Notably, there has been a substantial reduction in the values of RMSE, MAE, and MSE, while R ² has increased from 0.65 to 0.93. Figure 4C illustrates the comparison between actual and predicted values.

The performance evaluation of CART tree regression models, including RMSE, MAE, R ², MSE and training time T, is presented in Table 4. The iterative results of the hyperparameter optimization algorithm are illustrated in Figure 5B. It can be observed from the figure that the algorithm achieves convergence after three iterations with optimized hyperparameters: a minimum leaf size of 2.

The optimization of the model has led to a noticeable improvement in performance, as evident from the data presented in Table 4. Specifically, there has been a reduction in RMSE, MAE, and MSE values, indicating enhanced accuracy. Moreover, the R ² value has increased significantly from 0.84 to 0.89. Figure 4D illustrates the correlation between actual and predicted values.

The performance metrics, including RMSE, MAE, R ², MSE, and training time T of the boosting tree and random forest regression models in the ensemble tree regression approach are presented in Table 4. The iterative results of the hyperparameter optimization algorithm are illustrated in Figure 5C. It can be observed from the figure that the algorithm converges at the 23rd iteration with optimized hyperparameters as follows: Boosting algorithm is employed for ensemble learning method with a total of 497 learners, minimum leaf size set to 7, learning rate set to 0.0675, and number of sampled predictors limited to 5.

The performance of the empirically defined random forest regression model is superior to that of the boosting tree regression model, as evident from Table 4 prior to hyperparameter optimization. Moreover, the random forest regression model exhibits smaller values for RMSE, MAE, and MSE compared to the boosting tree regression model. Additionally, the R ² value of the random forest regression model surpasses that of the boosting tree regression model by 0.01. Subsequent optimization resulted in decreased RMSE, MAE, and MSE values along with an increased R ² value of 0.94, leading to a more refined and optimized model. Figure 4E illustrates the comparison between real and predicted values.