Predicting the sheet resistance of laser-induced graphitic carbon using machine learning

While laser-induced graphitic carbon (LIGC) has been used to fabricate cost-effective conductive carbon on flexible substrates for applications such as sensors and energy storage devices, predicting the resistance of the component fabricated via LIGC remains challenging. In this study, a two-step machine learning-based modeling framework is developed to predict the sheet resistance of the materials fabricated using LIGC. The two-step modeling framework consists of classification and regression. First, random forest (RF) is used to classify successful and failed trials. Second, extreme gradient boosting (XGBoost), RF, support vector machine with radial basis function, multivariate adaptive spline regression, and multilayer perceptron are used to predict the sheet resistance in each successful trial. In addition, an analysis of the change in sheet resistance with respect to laser energy per unit area is conducted to remove data points with high sheet resistance. XGBoost is also used to determine the importance of each process parameter. We demonstrate the modeling framework on datasets collected from experiments where LIGC lines (1D) and LIGC squares (2D) are engraved. For the 1D dataset, the RF classification model achieves a 95% accuracy. For both 1D and 2D datasets, a comparative study shows that XGBoost outperforms other algorithms. XGBoost predicts the sheet resistance of the LIGC lines and squares with a MAPE of 7.08% and 8.75%, respectively. XGBoost also identifies laser resolution as the most significant parameter. Moreover, experimental results show that models built on the dataset merging the 1D and 2D datasets result in lower prediction accuracy than those built on the 1D and 2D datasets separately. The modeling framework allows one to determine the sheet resistance of LIGC with varying laser processing conditions without conducting expensive and time-consuming experiments.


Introduction
The fabrication of graphene is crucial to a wide range of applications, including biomedical sensors, energy storage, and electrocatalysis due to its high thermal and electrical conductivity, strong mechanical properties, and biocompatibility [1].The scalable deposition and patterning of graphene over large areas on flexible substrates can enable wearable and flexible electronics systems [2,3].However, traditional top-down or bottom-up graphene fabrication techniques require harmful chemicals or hightemperature vacuum processes, which are complex, expensive, and potentially incompatible with flexible substrates [4].Furthermore, the resulting material needs to be patterned by photolithography and etching, which are also complex, multi-step processes that require the preparation of photomasks [5].By contrast, graphene can be fabricated by directly converting polymeric substrates into graphene using a laser [6,7].This one-step process, known as laser-induced graphene (LIG) or laser-induced graphitic carbon (LIGC), does not require separate patterning or post-processing.LIGC is compatible with flexible substrates, most commonly polyimide, 3D printed polymers [8], and other organic materials, including paper and wood [9].The resulting material has a porous morphology with a large specific surface area, making it ideal as electrodes for micro-supercapacitors.Depending on the laser conditions and polymer precursor, the material properties can vary, including flake size, number of graphene layers in each flake, and concentration of non-carbon atoms [10,11].Different devices and applications require different LIGC properties; however, the sheet resistance of LIGC is a fundamental property that affects the performance of most electrical LIGC-based devices in different applications.It is, therefore, essential to understand how laser conditions affect the sheet resistance.
Many studies have been conducted to optimize the LIGC sheet resistance based on the underlying physics of the process.Minhas-Khan et al [12] discussed several studies on polyimide with sheet resistance in the range of 10-50 Ω sq −1 .By increasing the overlap between adjacent pulses and maximizing the energy delivered to the polymer, a sheet resistance of 6.14 Ω sq −1 was achieved with a CO 2 laser on polyimide.Tavakkoli Gilavan et al [8] created very thick LIG on 3D-printed polyetherimide, which resulted in a very low sheet resistance of 0.30 Ω sq −1 .Morosawa et al [13] proposed a LIG process with the femtosecond laser and cellulose nanofiber film substrate to fabricate highly conductive graphitic carbon with a conductivity of 6.9 S cm −1 .For biomedical applications, Zhu et al [14] developed a coating process of various uniform metals on the LIG sheet with electroless plating for sensing glucose in diabetes patients.A good conductance of 25 S sq −1 was achieved.Joanni et al [15] introduced a novel method of transferring the graphene pattern through glass plates to arbitrary substrates.With this process, a LIG pattern was successfully transferred from polyimide to polydimethylsiloxane, and a low resistivity of 68 µΩ m −1 was achieved.Nova and Zarzar [16] introduced a direct laser writing method from liquid organic precursors.An inorganic glass substrate, CO 2 laser, and liquid precursors were utilized to fabricate sheets with sheet resistance as low as 4 Ω sq −1 .
In addition to these experimental studies, the LIGC process was investigated using numerical simulations.While multiphysics models were developed to understand the temperature evolution during the LIGC process [12], further investigations are needed to give direct insights into the conversion of the polymer into LIGC.Molecular dynamics simulations were used to study the LIGC process at the atomic level [17,18]; however, this approach was limited in terms of the spatial and temporal scales that could be simulated due to computational limitations.Machine learning is an alternative approach to modeling the LIGC process.Using Bayesian model-based optimization, a notable improvement of the Raman G/D ratio, which indicates the degree of graphitization in the LIG patterns, was achieved by a factor of four [19,20].Furthermore, to monitor the process of LIG formation, computer vision and deep transfer learning models were developed [21].Unlabeled images of LIG sheets were classified with labeled images as training data; thus, graphitization quality could be evaluated during manufacturing.An artificial neural network model was developed to track the degradation and predict the life cycle of LIG integrated into a composite structure.The non-linear damage progression within the structure was realized [22,23].Despite the advancements in applying monitoring and predictive models, predicting LIG sheet resistance using machine learning has yet to be reported.
To fill this research gap, this paper aims to develop a two-step machine learning-based framework to predict the LIGC sheet resistance.
The main advantage of machine learning models is the ability to discover hidden, complex, and non-linear relationships between independent variables and target variable(s) when the knowledge of the underlying physics of a process or a system is limited.The two-step machine learning-based framework we developed can predict the LIGC sheet resistance with different laser configurations and determine the importance of process parameters.Due to the presence of non-numerical sheet resistance values when the polymer was not converted to LIG, both classification and regression algorithms were implemented.First, a machine learning algorithm was used to classify these non-numerical values indicating success or failure in producing graphitic carbon on the material surface due to insufficient or excessive laser power.Second, another machine learning algorithm was used to predict the LIGC sheet resistance for successful engravement.Random forest (RF) was chosen to build the classification model.Extreme gradient boosting (XGBoost), RF, support vector machine with radial basis function (SVM-RBF), multivariate adaptive spline regression (MARS), and multilayer perceptron (MLP) were used to build regression models.These algorithms were selected due to their excellent performance, robustness, and wide applications in various manufacturing fields [24][25][26][27][28].However, due to outliers and noise in the raw data, we developed a filtering method by analyzing sheet resistance change with respect to laser power.The contributions of this work are as follows: • This study demonstrated that machine learning models can predict the sheet resistance of the LIGC process with high accuracy.• Important laser parameters, which significantly affect the sheet resistance, were identified using XGBoost.

Experimental setup and data description
The fabrication and characterization of the LIGC sheet resistance were comprehensively described and analyzed by Minhas-Khan et al [12]; hence, only a summary of the LIGC sheet resistance data collection is provided.This section focuses on the discussion of the independent variables used to build the models.In short, a 75 W CO 2 laser (Epilog Fusion M2 Laser) was used to convert a 1 cm by 1 cm square area of commercially available polyimide tape (Kapton, 125 µm thickness) into LIGC.The laser power was varied by pulse width modulation as a percentage of the maximum 75 W.The sheet resistance of the LIGC was measured using the four-point resistance measurement method.Figure 1 shows the experiment procedures for engraving different LIGC configurations.Two datasets were generated from the experiments.One is referred to as 1D (LIGC lines) dataset where LIGC lines were engraved; the other is referred to as 2D (LIGC squares) dataset where LIGC squares were engraved.The 1D and 2D datasets contain 200 and 370 data points, respectively.The 1D dataset contains information about (1) laser lines (number of parallel lines that are engraved or irradiated), (2) laser resolution (dot per inch), (3) laser velocity (millimeter per second), (4) a percentage of 75 W laser power (%), (5) energy per pulse (joules), (6) energy per unit area (EPUA) (joules per centimeter squared), (7) irradiated width (micrometers), and (8) sheet resistance (Ω sq −1 ).The 2D dataset contains similar information without laser lines and irradiated width.
In the 2D experiments, a square (or grid) pattern was used to guide the laser motion; hence, the number of laser lines passing over any point at the center of the square is directly proportional to laser resolution, and irradiated width is constant and practically infinite (4 mm by 4 mm).While there are eight and six variables in the 1D and 2D datasets, respectively, laser EPUA and laser energy per pulse were not considered when building predictive models because they were determined based on other variables.Another characteristic of the dataset was the non-numerical values of sheet resistance (i.e., phrases or string values) when no conductive LIGC was created, which were recorded as 'ablated' or 'no LIGC engraved' .These nonnumerical values, which were observed only in the 1D dataset, describe two thresholds for laser power.Laser power exceeding the upper threshold ablated the material; hence, 'ablated' data points were recorded.Insufficient laser power, which was below the lower threshold, could not produce graphitic carbon.'no LIGC engraved' indicates this phenomenon as engravement was not produced on the material surface.From a physics perspective, these data points describe different experimental outcomes, and they can all be labeled as 'failed' (i.e., failure to produce any engravement) to build classification models.In the 1D experiments, laser scanning velocity was constant, whereas the number of laser lines, laser resolution, irradiated width, and laser power varied.In contrast, in the 2D experiments, the irradiated width was constant, whereas laser resolution, velocity, and power varied.The number of laser lines is directly proportional to laser resolution, which is not an independent parameter for the 2D dataset.Tables 1 and 2 provide more details about the experimental designs for the 1D and 2D datasets.Data distributions for the 1D and 2D datasets are shown in figures 2 and 3.

Predictive modeling framework
Figure 4 illustrates the overall predictive modeling framework.The evaluation metrics include accuracy for the classification model and mean absolute percentage error (MAPE) for regression models.After being optimized with a tuning algorithm on the training dataset, the hyperparameters were tested on the test dataset, and the MAPE on the test datasets was calculated.The machine learning algorithms were tested on three train/test split ratios, including 70/30, 80/20, and 90/10.

Hyperparameter tuning
Tuning hyperparameters is crucial for building accurate machine learning models.The two most commonly used tuning techniques are grid search and random search.The grid search method divides hyperparameters into a grid space.Each point within the grid space represents a combination of the hyperparameters.The grid search method is a brute force method; hence, its computational expense scales with the size of the dataset.The random search method operates similarly to the grid search method; however, instead of searching all possible combinations in the grid space, random search only tests some combinations randomly.Thus, a random search may only find suboptimal hyperparameters.In this study, Optuna was used to address the limitations of traditional tuning methods.Optuna is an effective hyperparameter tuning method implemented to tune various machine learning algorithms [29].The backbone of Optuna is the tree-structured Parzen estimator algorithm (TPE) [30].To find the optimal hyperparameters, TPE uses Bayes rule, and it consists of a series of steps: (1) define the range of hyperparameters and record some initial results with random combinations of hyperparameters, (2) form two groups from the initial records based on their quantile, then calculate their densities, (3) model the two-group based on Parzen estimator to suggest a probability of promising hyperparameters that likely yield optimal results.TPE is an iterative process; subsequent steps are suggested based on the density functions of previous results.Moreover, Optuna also features pruning (or early stopping) which allows the tuning algorithm to stop searching for hyperparameters in specific directions if the results do not improve within specific iterations.Overall, compared to the grid and random search, Optuna is a more effective tuning algorithm for this study due to its informed searching sequences.

Building a classification model with random forest
RF was selected to perform the classification task due to its unique advantages.Tree-based methods such as RF are not sensitive to data scaling and multicollinearity; hence, cumbersome data preprocessing can be avoided.RF randomly selects subsets of features to create trees; then, all the trees are bagged.Each tree represents a vote, and the initial classification of an input vector is made based on the most votes.The decisions on choosing the nodes are based on information gain.RF makes decisions to maximize information gain.In this work, entropy was used as information gain, which is defined as: and entropy (or log loss) is defined as: where x i is the instance (observation), in feature (row) A, p i indicates the probabilities of x i being classified in C i (where C i is the correct label of x i ).Info (x i ) is the entropy of the parent node (before splitting), and Info A (x i ) is the entropy of the child node (after splitting).The tree algorithm implemented in RF is classification and regression tree which only creates a binary tree to find the best-split node for one independent variable at a time [31].This tree algorithm may lead to an unstable tree causing poor generalized prediction performance.Thus, hyperparameter tuning must be carefully executed to combat this overfitting issue.RF models were built to classify successful and failed trials indicating the conditions of LIGC engravement.Stratified sampling was implemented to split datasets into three train/test split ratios, including 70/30, 80/10, and 90/10.This concludes the explanation and model-building steps for the RF classification models.

Preprocessing data for regression models
Data preprocessing consists of several steps.First, we implemented standardization.Second, we created a   new dataset by combining 1D and 2D datasets.Third, we removed data points with high resistance values based on their changes with respect to laser EPUA.Standardization is a scaling method where the values are centered around the mean with a unit standard deviation.While standardization does not affect XGBoost and RF models, MARS, SVM-RBF, and MLP could benefit from standardization since they are affected by scaling.For instance, the MARS model is an ensemble of multiple piecewise linear functions characterized by hinge functions.Linear functions describe the relationship between independent and dependent variables with coefficients.As such, largescaling variables have more influence than smallerscaled ones; however, this does not necessarily mean they are significant predictors.As a result, standardization was performed on independent variables for the 1D and 2D datasets to ensure the consistency of data scales, which theoretically should enhance the performance of the MARS, SVM-RBF, and MLP models.
A new dataset was created by combining 1D and 2D datasets because they contain LIGC resistance data collected from similar experiments where lines and squares were created.To merge the 1D and 2D datasets, the irradiated width was removed from the 1D dataset.The number of laser lines was added in the 2D dataset to signify that the number of overlapping laser passes is significantly larger in squares than in lines.However, in large squares, the total number of lines is directly proportional to DPI and not an independent parameter.Furthermore, it does not have the same physical meaning as in 1D lines since the number of laser passes in every position does not depend on it.Therefore, an arbitrary dummy value can be chosen for the number of lines in large squares.Based on the experimental data, the number of lines for the 2D dataset ranges from 24 to 189.We used 189 as the dummy variable.
We also filtered high resistance values based on their changes with respect to laser EPUA.The initial MAPE for all regression models was relatively high without filtering out these high resistance values.The high MAPE was mainly because the models could not predict high resistance values.Examining the plots of the LIGC sheet resistance with respect to laser energy per pulse (EPP) and laser EPUA in figure 5 reveals a pattern: the resistance for each experimental configuration almost always starts at a significantly high value before stabilizing to a semi-linearly path.Since this study aims to predict the low resistance value achieved by maximizing energy delivery and pulse overlap, these high resistance values are generally not of primary interest; thus, the data points with high resistance values are considered outliers.The EPUA and energy per pulse equations are described as follows: where power denotes the percentage power of the 75 W CO 2 laser, resolution denotes the laser resolution value in dots per inch, and velocity denotes the laser velocity value in millimeters per second.This study proposes an analysis of change in resistance with respect to EPUA, which is referred to as resistance gradient for short.The resistance gradient was calculated by dividing the change of sheet resistance by the change of EPUA.Intrinsically, sorting the resistance values in an ascending order will assist in removing these high resistance values.However, since sheet resistance data points were recorded within their specific experimental designs, each resistance value should be compared to others within its experimental groups rather than the entire population.The datasets were divided into different groups to determine the resistance gradient.For the 1D dataset, three groups were created based on the number of laser lines.Each group was further divided into three subgroups based on the laser resolution.By contrast, the 2D dataset was divided based on different criteria.Since the 2D dataset contains only square laser patterns, the first criterion was laser resolution, which resulted in four groups.However, while sorting and evaluating the resistance gradient for the 2D dataset, errors of dividing by zero change in total laser EPUA occur, i.e., two or multiple data points with equal EPUA.For instance, the EPUA value of a data point with a laser velocity of 15.875 mm s −1 and a percentage laser power value of 14% is identical to one of a data point with a laser velocity value of 3.175 mm s −1 and a percentage laser power value of 7%.Therefore, even though some values in the 2D dataset share identical values of EPUA, they describe different experimental parameters.Hence, each laser resolution group in the 2D dataset was further divided based on the laser velocity.A summary of the group formation process is shown in table 3.
In summary, the 1D dataset was divided into 12 groups of four main groups (laser lines) and three subgroups (laser resolutions).The 2D dataset was divided into 16 groups of four main groups (laser resolutions) and four main subgroups (laser velocity).Figure 5 shows that different experiment groups exhibited varying ranges of resistance.Since characterizing the resistance gradient to filter high data points requires the resistance gradient values between all experiment setups to be comparable, a normalization step, which divides each calculated resistance gradient by the difference between the highest resistance and lowest resistance within each group, was applied.Figures 6 and 7 show the original values of the normalized resistance gradient in the 1D and 2D datasets.The plotted results are scaled differently based on the range of the groups.
Based on figures 6 and 7, high resistance gradients in the 1D and 2D datasets were observed.For the 1D dataset, the highest resistance gradient was close to −40 (i.e., the normalized resistance decreases at 40 per J cm −2 ), followed by a value of −28.For the 2D dataset, the highest resistance gradient of more than −250 was observed.Removing these exceptionally high resistance data points can improve the performance of the predictive models.The values of ±10 were used as the thresholds to remove outliners for the 1D and 2D datasets, respectively.Figures 8 and 9 show that the maximum resistance gradients were approximately −6 and −10 for the 1D and 2D datasets, respectively, after removing the high resistance data points.Noticeably, some lower resistance data points were also removed.While these resistance values are low compared to the entire population, they are exceptionally high resistance data points within their groups.Hence, these data points can only be detected filtered through the gradient process.
In addition, the 1D and 2D datasets are recognized as small datasets, which may be subjected to disproportionate splitting and overfitting.Disproportionate splitting occurs when the distribution of samples in the training data and test data does not represent the population distribution.A disproportionate sampling may produce a low bias accuracy.For instance, a random sampling process may output a test dataset consisting of a large fraction of low-value points and a negligible sample for middle to high-value points.High-accuracy regression models built to predict this test dataset can only predict low-value data points.These models are impractical due to the distortion of the sampling training and test dataset.Hence, in this study, stratified sampling was used to address this issue.Moreover, overfitting is another issue when building predictive models with small datasets.The smaller the dataset, the more models can be built to fit.Complex models may yield excellent accuracy on training data while returning a poor performance on test data.This study performed hyperparameter tuning with cross validations on the dataset to prevent overfitting.

Building a regression model with XGBoost
XGBoost is a tree-based algorithm utilizing gradient descent to minimize the loss function in an iterative sequence [32].XGBoost can be implemented for both regression and classification tasks; however, for this study, XGBoost was implemented only to build regression models.XGBoost is built upon the foundation of the traditional gradient-boosting algorithm proposed by Friedman [33].Like Friedman's gradient boosting algorithm, XGBoost still uses weak learners (i.e., decision stumps or trees with singular nodes) to build predictive models.However, compared to the vanilla version, XGBoost improves computing time by introducing second-order Taylor's series approximation for loss function and parallel processing.Moreover, regularization is also implemented in XGBoost to prevent the model from becoming overcomplex, leading to overfitting.Prediction outputs for XGBoost regression are defined as: where ŷ denotes the predicted LIGC sheet resistance values, x i are the input variables which are defined as input variables in 1D and 2D datasets (tables 1 and 2), f k = ω q (x) is the additive function evaluated as tree structure q and leaf weights ω.Equation ( 5) explains that the prediction output is calculated as the sum of the weight function in each weak learner.Like most other classification and regression algorithms in machine learning, XGBoost optimizes an objective function by minimizing a loss function to train   the model additively.The objective function used in additive training of XGBoost can be described as: where 1 is the loss function, which must be differentiable and convex for gradient descent, Ω is the regularization term, and t denotes the t-iteration.
Instead of the traditional optimization by searching the Euclidean space, the function describes the optimization process by adding a tree represented by the term f t (x i ).In the other term, the tree, which improves the objective function the most, will be added in the next iteration.The penalty term Ω is defined as: where T stands for the terminal, denoting the nodes (or leaves) in a tree, γ is the user-defined parameter used to prune the tree, thus, preventing over-complexity, and λ (lambda) is the penalty term that shrinks the leaf weights.Like ridge regression, if λ = 0, regularization is not being implemented; if λ > 0, leaf weight shrinks toward zero.To evaluate the objective function f t (x i ) term in equation ( 7), secondorder Taylor's series approximation is used, and the objective function can be re-written as follows: where ) and ) are the first and second derivatives of the loss function f t (where g i stands for gradient and h i stands for Hessian).Removing the constant terms (i.e., terms do not contain leaf weight ω).The objective function can be re-written as: to find the optimal value for L (t) , one can take the first-order derivative of the objective function and then set it to zero.Solving for ω yields: As mentioned, XGBoost optimizes the objective function by adding the optimal leaf weight score ω of the tree structure q.Then, scores measuring the quality of tree structure q can be evaluated to find the following optimal tree structure q.Define I j as the instance set of leaf j, this scoring metric of quality for a given tree q is defined as: In this study, the loss function for XGBoost was the squared error loss function.The squared loss function for a single observation and prediction is defined as: substituting f t as the squared loss function error and evaluating h i and g i , weight ω from equation (10) will be equal to Sum of residuals Number of residuals+λ .Hence, XGBoost regression with squared error loss function trains the model iteratively by minimizing the objective function through preceding residuals.Since this study built predictive models on relatively small datasets, XGBoost must be implemented cautiously to avoid overfitting.XGBoost employs exact greedy algorithms to search for its optimal learner; that is, the algorithm picks the best possible split loss reduction after iterating over all given input variables.This behavior leads to the issue where the algorithm choice for local optimality might not equate to an acceptable global solution.Thus, tree complexityrelated and sampling parameters, such as tree depth, tree pruning, regularization parameters, etc., were carefully tuned to ensure the model is conservative enough to combat overfitting.As a reminder, five-fold cross-validation was also used to combat overfitting of the predictive model.The results and discussion section will provide further information for XGBoost hyperparameters.

Results and discussion
The 1D dataset for classification contains 89 and 111 'success' and 'fail' data points.Figure 10 shows the classification accuracy of RF.As expected, reducing the size of the training data size reduced prediction accuracy.
XGBoost models were built for 1D and 2D datasets.Figure 11 shows the prediction accuracy of XGBoost for the 1D dataset in terms of MAPE for train/test split ratios of 90/10, 80/20, and 70/30.The best prediction accuracy of 7.08% was achieved by XGBoost on the 90/10 split ratio.The predictive performances of XGBoost decreased as the training dataset size decreased: the model trained with the 90/10 split ratio achieved the best accuracy, followed by the one trained on the 80/20 split ratio as the second best, and the models trained on the 70/30 split ratio as the least accurate one.Predictions for low LIGC sheet resistances (lower than 100 Ω sq −1 ) were relatively accurate.However, XGBoost still struggled to capture the high LIGC sheet resistance data points.This phenomenon could be attributed to the disproportionate distribution of the data points for the low and high-value LIGC sheet resistance in the 1D dataset: the population of low-value data points is much denser than the high-value data points.Compared to the prediction accuracy for high resistance data points, low-value data points would likely to be predicted with higher accuracy since more information about the low-value data points is readily available to learn.If these high sheet resistance values are considered as outliers recorded while increasing delivered EPUA, they can be filtered out during the resistance gradient filtering process by adjusting the thresholds.In this case, the performance of XGBoost algorithms should increase significantly since the high MAPE was mainly caused by the prediction errors of these high resistance values.
Figure 12 shows the prediction accuracy of XGBoost for the 2D dataset.The best model was achieved with the 90/10 split ratio (8.75% MAPE).Similar to the 1D dataset results trend, the accuracy was lower when the training dataset size reduced.Another similar trend can be observed in 2D dataset results: low resistance values (lower than 100 Ω sq −1 ) were predicted with high accuracy.In comparison, the prediction errors for high resistance values (higher than 100 Ω sq −1 ) were much higher.Thus, improved accuracy can be achieved by adjusting the filtering thresholds in the resistance gradient analysis, which could remove these high resistance values.
Figure 13 shows the prediction accuracy of XGBoost for the merged dataset.The XGBoost model built with the 90/10 split ratio consistently achieved the lowest MAPE at 14.63%.Figure 13 suggests that the prediction accuracy of XGBoost on the merged dataset is worse than those on 1D and 2D Tables 4 and 5 show comprehensive comparisons.Tables 6 and 7 provide the hyperparameters for the RF classification and XGBoost regression models with optimal accuracy and MAPE.
Figure 14 shows a summary of all XGBoost regression models and the feature importance identified by XGBoost models.XGBoost achieved excellent prediction accuracy on the 1D and 2D datasets.In       14 shows that for 1D and 2D datasets, laser resolution was the most significant parameter when building XGBoost predictive models.The finding that laser resolution was the most significant feature is consistent with experimental and numerical simulation results.These results have demonstrated that the quality of LIGC increases with larger DPI values due to the more closely spaced pulses, smaller energy per pulse, and a larger number of consecutive pulses for every point on the surface, which leads to a more consistent temperature distribution [12].Table 5 suggests that while the MARS models predicted LIGC sheet resistance on the 1D dataset with relatively high accuracy, the MARS models on the 2D dataset could potentially be under-fitted due to the large MAPE; thus, increasing models' complexity could theoretically increase models' performance on the test dataset.However, for the MARS models trained on the 2D dataset, MAPE for five-fold cross validation results on the training datasets were approximately equal to the ones evaluated on the test datasets.This result suggests that the models were not underfitted, but failed to learn the trend of the 2D training and test datasets.Perhaps, the MARS models suffered from multicollinearity [34].During the forward selection process, multicollinearity (or collinearity if MARS only considers two independent variables) forces the algorithm to place its knots on one of the independent variables.This choice is arbitrary if the loss function is roughly identical; however, this arbitrary decision dramatically affects the subsequent knots.The MAPE of the RF and SVM-RBF regression models on the 1D and 2D datasets were comparable to those of the MARS models.For the RF regression models, the results suggest that the models could not accurately predict the LIGC sheet resistance.While XGBoost implemented boosting and bagging to build models, the RF regression models were built through bagging only.Hence, it implies that bagging alone could not generalize the LIGC data.The SVM-RBF models also exhibited comparatively high MAPE on both datasets.These results may be attributed to the assumption of the Gaussian distribution in the SVM-RBF algorithm.MLP generally requires large volumes of data to build accurate models [35].Data manipulation techniques can improve MLP's accuracy [36][37][38].The performance of MLP models were relatively poor because the 1D and 2D datasets are relatively small.It should be noted that while the performance of traditional machine learning models (i.e., SVM-RBF, RF, MARS, and XGBoost) trained on the 1D dataset was superior to the performance of the models trained on the 2D dataset, the reversed result was observed for MLP.Another observation is that the performance of the XGBoost models trained on the merged dataset were worse than ones of the XGBoost models trained on the 1D and 2D datasets separately.This is perhaps because the data distributions for 1D and 2D datasets are different.

Conclusion
A two-step machine learning-based computational framework was introduced to predict the LIGC sheet resistance.For the 1D dataset containing data collected from a single line, two lines, and three lines of laser engravement, RF was implemented to build classification models to classify successfully engraved trials and failed ones.An accuracy of 95% was achieved with the classification model.Then, the comparative study showed that XGBoost outperformed MARS, RF, SVM-RBF, and MLP.The XGBoost models predicted the LIGC sheet resistance of the 1D and 2D datasets with a MAPE of 7.08% and 8.75%, respectively.Furthermore, XGBoost identified laser resolution as the most significant parameter, which is consistent with the experimental results reported in our previous study where we found that increasing laser resolution reduced sheet resistance.We also found that the prediction accuracy of the models built on the merged dataset was lower than the prediction accuracy of the models built on the 1D and 2D datasets separately.

Figure 1 .
Figure 1.Schematics of experimental procedures describing (a) different laser lines, (b) (i) single line of LIGC at different laser resolutions; (ii) different number of laser lines at fixed laser resolution; (iii) square LIGC pattern.Reprinted from [12], Copyright (2021), with permission from Elsevier.

Figure 5 .
Figure 5. Sheet resistance distribution for 1D dataset with respect to (a) energy per pulse and (b) energy per unit area.Representatives of sheet resistance distribution for the 2D dataset with respect to (c) energy per pulse and (d) energy per unit area.For each laser resolution group in the 2D dataset, three representatives of the velocity groups are plotted.Reprinted from[12], Copyright (2021), with permission from Elsevier.

Figure 8 .
Figure 8. Normalized resistance gradient for the three groups in the 1D dataset after applying filters for (a) one line group (48 data points), (b) two lines group (21 data points), and (c) three lines group (17 data points).

Figure 10 .
Figure 10.Confusion matrix for the RF classification model for three different train/test split ratios (a) 95% accuracy with 90% training dataset, (b) 90% accuracy with 80% training dataset, and (c) 83.33% accuracy with 70% training dataset.

Table 1 .
Design of experiments for the 1D dataset.

Table 2 .
Design of experiments for the 2D dataset.

Table 3 .
Summary of criteria for group formation.
Figure 6.Original normalized resistance gradient for the three groups in a 1D dataset with thresholds as red dash lines for (a) line group (51 data points), (b) two lines group (21 data points), and (c) three lines group (17 data points).

Table 4 .
Results for the comparative study of machine learning models for the 1D dataset in terms of MAPE (%).

Table 5 .
Results for the comparative study of machine learning models for the 2D dataset in terms of MAPE (%).

Table 6 .
Hyperparameters for the RF classification model with the train/test split ratio of 90/10.

Table 7 .
Hyperparameters for XGBoost with a train/test split ratio of 90/10.