Investigation of Melt Pool Geometry Control in Additive Manufacturing Using Hybrid Modeling

: Metal additive manufacturing (AM) works on the principle of consolidating feedstock material in layers towards the fabrication of complex objects through localized melting and resolidiﬁcation using high-power energy sources. Powder bed fusion and directed energy deposition are two widespread metal AM processes that are currently in use. During layer-by-layer fabrication, as the components continue to gain thermal energy, the melt pool geometry undergoes substantial changes if the process parameters are not appropriately adjusted on-the-ﬂy. Although control of melt pool geometry via feedback or feedforward methods is a possibility, the time needed for changes in process parameters to translate into adjustments in melt pool geometry is of critical concern. A second option is to implement multi-physics simulation models that can provide estimates of temporal process parameter evolution. However, such models are computationally near intractable when they are coupled with an optimization framework for ﬁnding process parameters that can retain the desired melt pool geometry as a function of time. To address these challenges, a hybrid framework involving machine learning-assisted process modeling and optimization for controlling the melt pool geometry during the build process is developed and validated using experimental observations. A widely used 3D analytical model capable of predicting the thermal distribution in a moving melt pool is implemented and, thereafter, a nonparametric Bayesian, namely, Gaussian Process (GP), model is used for the prediction of time-dependent melt pool geometry (e.g., dimensions) at different values of the process parameters with excellent accuracy along with uncertainty quantiﬁcation at the prediction points. Finally, a surrogate-assisted statistical learning and optimization architecture involving GP-based modeling and Bayesian Optimization (BO) is employed for predicting the optimal set of process parameters as the scan progresses to keep the melt pool dimensions at desired values. The results demonstrate that a model-based optimization can be signiﬁcantly accelerated using tools of machine learning in a data-driven setting and reliable a priori estimates of process parameter evolution can be generated to obtain desired melt pool dimensions for the entire build process. This paper develops a novel hybrid methodology for control of melt pool geometry in AM in the framework of ML-assisted modeling and optimization. The continuous changes in the melt pool geometry are predicted by a low-cost GP surrogate developed using an experimentally validated analytical 3D model. The uncertainties in the predictions for the melt pool geometry are quantiﬁed using GPs. Reliable estimates of the optimal process parameter evolution are obtained through active learning via BO by devising appropriate objective functions that quantify the control requirements. The methodology provides an estimation of process parameter variation during the AM deposition process in order to maintain a desired melt pool geometry under continuously varying thermal conditions.


Introduction
Metal additive manufacturing (AM) facilitates direct fabrication of near-net-shape metallic components, prototypes, or both under rapid solidification conditions [1]. AM is currently used in manufacturing a wide variety of components with increasing complexity, for example, fuel nozzles, rocket injectors, and lattice structures [2]. The concept of AM is built on the principle of incremental layer-by-layer material consolidation through localized melting and resolidification of feedstock materials by using high-power energy sources [3]. The localized heating causes the formation of a melt pool that controls the microstructure and, therefore, the properties of the manufactured component [3]. Due to the cyclic nature of the deposition process as AM components continue to gain thermal energy, the thermal gradient (G) and the solidification velocity (R) inside the melt pool continuously change, resulting in significant alterations in the melt pool properties (e.g., geometry, thermal profile, and flow field among others) between the initial and final layers [4,5] as illustrated in Figure 1a. Figure 1b shows that in spite of fixing the G/R ratio at a value meant for yielding a columnar microstructure, at the leading edge (e.g., beginning of the scan), the grains are columnar; however, at the trailing edge (e.g., towards the end of the scan), the grains become equiaxed [6] due to progressive heating. The primary reason for controlling melt pool geometry in metal AM is to allow a part to be built with consistent melt pool dimensions, even as thermal conditions change continuously [7]. Understanding the fundamental physics of the melt pool evolution is, therefore, a key requirement for AM process development, optimization, and control. Although high-fidelity computational modeling of AM processes can provide reliable estimates of the melt pool evolution during the build process, such elaborate models are nearly-intractable when coupled with an optimization code for melt pool geometry control due to the computational costs [8]. Therefore, while these models are well suited for understanding the physical phenomena, challenges exist in using them for performing process design, optimization, and control [9].
In this respect, machine learning (ML) shows potential in assisting and automating the process of manufacturing [10] through reliable prediction of melt pool geometries [11]. For example, Neural Networks (NN) [12] are widely used as a popular choice in prediction problems by modeling nonparametric input-output relationships. Despite the simplicity of usage, most of these techniques based on complex NN architectures suffer from issues of interpretability [13]. Moreover, the applications of such methods in limited datasets are rare, and lack of data can often result in poor predictive models that lack generalizability due to overfitting [14]. Moreover, a vast majority of the ML techniques used in AM relies on point estimates of the quantities of interest, without catering for uncertainty in the predictions. In critical applications involving high stakes associated with mispredictions, the estimates of uncertainty are particularly important. An attractive alternative, therefore, is to use probabilistic ML techniques like Gaussian Processes (GPs) [15] that offer the advantages of interpretability and applicability in limited data regimes. However, their applications in surrogate-based modeling and optimization are relatively less explored in multi-physics problems in AM. An implementation of surrogate modeling through the construction of computationally efficient approximations that can be used in lieu of the original simulation model [16], therefore, holds significant promise in metal AM.  [4,5], with permissions of Elsevier, 2009 and Springer Nature, 2016). (b) Simulation results showing transition from columnar grains at the beginning of the scan to equiaxed grains towards the end of the scan in an Inconel 718 specimen with a G/R ratio meant to yield columnar microstructure (Reproduced from [6], with permission of Elsevier, 2019). Surrogate-assisted modeling and optimization techniques are popularly used in various applications that involve expensive computational models [17][18][19][20][21][22][23] for function evaluations. Among the different types of surrogate models, GPs [15] are used profusely for modeling black-box functions whereby a fully Bayesian approach allows for probabilistic estimates of the target functions [24][25][26][27].
With a relatively small number of measurements, a GP surrogate can be learnt to serve as a proxy to an expensive objective function [28] (e.g., prediction of melt pool dimensions for a range of process parameters in metal AM [29]). Under the settings of a GP surrogate, a Bayesian Optimization (BO) set-up can be invoked for gradient-free global optimization of an objective function under budget limitations [30] (e.g., prediction of optimal process parameters for controlling temporal melt pool dimensions in metal AM within N number of iterations; N being an user input). Moreover, with nonlinearities in the objective function, a search for optimum would require significant amount of sampling in the search space, particularly in high dimensions. In such settings, BO is found to be quite successful [31].
In the field of AM, there are very few applications of GPs as surrogate models for expensive experiments and simulations. For example, Tapia et al. [32] used experimental data to learn spatial GPs for predicting porosity in metal AM produced during the laser powder bed fusion process. In another work, Tapia at al. [16] demonstrated the usage of GPs in predicting melt pool depth as a function of different process parameters, which were used to describe regimes of operation where the process was expected to be robust. Seede et al. [33] used a GP framework to develop a calibrated surrogate model for predicting optimal process parameters for building porosity free parts with low-alloy martensitic steel AF9628. However, to the best of the authors' knowledge, no work exists in the field of AM that leverages the surrogate-based predictions and uses them as a basis for active learning strategies for optimizing a desired objective, e.g., controlling melt pool dimensions over time under computational budget constraints in laser powder bed fusion process. The design and deployment of such predictive methodologies will immensely augment the response time of feedback or feedforward control strategies as the current response times are rather long compared to the process time scales [34].
To address this gap, this paper proposes a novel framework of controlling the melt pool geometry in laser powder bed fusion process by formulating it as an optimization problem through the integration of the tools of physics-based analytical modeling and data-driven analysis. The cardinal contribution of this work in AM is to bridge the gap in model-assisted prediction and control of melt pool geometry by using ML techniques that can potentially accelerate the process of melt pool geometry control. The results demonstrate that by using a low-cost surrogate-assisted modeling using GP and BO, it is possible to obtain an excellent estimation of process parameter evolution as a function of time in order to maintain the desired meltpool geometry (e.g., dimensions) throughout the build process. The framework consists of three critical steps: (i) evaluation of the thermal field using an experimentally validated 3D analytical melt pool evolution model which serves as a source of data for formulating a GP surrogate, (ii) development of GPs with flexible kernel structures capable of handling anisotropies in the dataset, and (iii) establishment of a BO framework involving the GP surrogate that aims to solve a global optimization problem in a gradient-free setting under a constrained budget. The computational budget is predefined in most practical problems. With an appropriate selection of initial design points and acquisition function for active learning of optimal design points, this work achieves estimates of process parameters with limited model iterations. The GP surrogates, being based on Bayesian inference, offer principled estimates of uncertainty in predictions. While the present work uses a 3D analytical model for demonstrating the efficacy of the proposed approach, the framework described can be applied to any other prediction models of users' choice.
The paper is organized in five sections including the current one. Section 2 describes the development of a hybrid modeling framework consisting of the 3D analytical model and its validation against experimental data. Section 3 presents the results and discusses the significant findings. The paper is summarized and concluded in Section 4 along with recommendations for future research.
Supplemental information in Appendix A provides the mathematical background for ML aspects of modeling and optimization.

Simulation Model
A well-known analytical model developed by Eagar and Tsai [35] that solves for the 3-dimensional temperature field produced by a traveling distributed heat source on a semi-infinite plate is used in the present work. This model has been previously used by several different researchers for evaluating melt pool evolution in laser powder bed fusion (L-PBF) [33] and directed energy deposition (DED) [36] AM processes when extensive evaluation of process parameter space is required. It is a distributed heat source modification of the Rosenthal's [37,38] solution for the temperature distribution produced by a traveling point heat source. As compared to Rosenthal's solution, Eagar-Tsai's model provides a significant improvement in prediction of temperature in the near heat source regions. Figure 2 explains the coordinate system used in the model. The heat source is traveling with a uniform speed of v in the X-direction, and is assumed to be a 2D surface Gaussian: Here, Q(x, y) is the power distribution per unit area on the X-Y plane of the specimen produced by the Gaussian beam of peak power P with a distribution parameter σ. According to Eagar-Tsai's model, the temperature T(x, y, z, t), at a particular location (x, y, z) and time t, with the initial temperature of the substrate being T 0 , is denoted as Here, α L is the absorptivity of the laser beam, a k ρc is the thermal diffusivity, ρ is the density, and c is the specific heat capacity of the material of the specimen. The primary assumptions of the model include the following.

1.
Convective and radiative heat transfer from the substrate to the environment are ignored.
As the process deals with metals which are good conductors, heat transfer through radiation is negligible [39] as compared to that due to conduction.

2.
The temperature dependence of the thermo-physical properties is not taken into account.

3.
The substrate is semi-infinite, therefore the increase in the surface temperature T 0 with time is negligible.

4.
Phase change of the material is not taken into consideration.
Eagar-Tsai's model, despite its limitations (e.g., its inability to model keyholing effect [33]), is widely used in literature for parameter space exploration through design of experiments (DoE)-based approaches. Given the wide range of process parameters in metal AM machines (e.g., L-PBF has over 130 parameters that can affect the final part quality [40]), such a DoE approach would require an exponential number of samples in the dimension of the parameter space to fully explore the design space. This makes it prohibitively expensive to explore the design space for building complex parts through experimentation or high-fidelity simulations. Computationally efficient surrogates can reduce the computational burden by a significant amount. Eagar-Tsai's model, when appropriately calibrated, provides one such alternative as a low-cost simulation model that are extensively used by several researchers in recent times [33,41,42]. This model especially works well for single-track and single-layer AM depositions. Owing to these characteristics, the validation and control experiments outlined in this paper are all single-track and single-layer.

Experimental Validation of the Simulation Model
Eagar-Tsai presented experimental validation of their model with a range of process parameters in welding carbon steel plates [35]. Other researchers validated this model for a range of different materials such as nickel- [43], iron- [33], and titanium-based [44] alloys among others to obtain satisfactory accuracy with the experimental data across a range of AM processes, e.g., L-PBF and DED [45], in recent times. However, most of the experimental validation of this model, so far, has been performed under steady state conditions. Being a transient model, Eagar-Tsai's formulation can potentially be used to simulate the melt pool dimensions as a function of time as well. In this section, both steady and unsteady state validation of the Eagar-Tsai's formulation are presented.

Steady-State Experimental Validation
For the steady-state model validation, the melt pool dimensions calculated by Eagar-Tsai's model are compared with those obtained using finite-element simulations and experiments for a popular nickel-based superalloy CMSX-4 reported by Wang et.al [36]. This alloy was processed by several different researchers using L-PBF [46,47], electron beam powder bed fusion (E-PBF) [48], and DED [49], and therefore a considerable amount of experimental data is  Table 1. A CO 2 laser beam with radius W was used in their work. Table 1. Processing parameters as reported in Wang et al. [36]. Figure 3 shows the comparison between the current work with those obtained experimentally and using finite element simulations by Wang et al. [36] showing excellent agreement. Figure 4 shows a transverse cross section (Y − Z plane, with x = 0) of the analytically calculated melt pool geometry juxtaposed on the fusion zone produced by the process parameters in the CMSX-4 Specimen A, as reported by Wang et al. [36]. Similar agreements are also observed for the other two specimens (i.e., B and C), as evident from the results reported in Figure 3. Cross-sectional images for specimens B and C are omitted for brevity.

Transient Experimental Validation
After obtaining excellent agreement with the steady-state results, the analytical model is further validated using transient melt pool results available in the open literature [51]. Figure 5 shows the longitudinal cross section of an L-PBF-processed CMSX-4 specimen. This specimen is fabricated by consolidating a powder layer thickness of 1.4 mm on a CMSX-4 substrate having dimensions of 35.56 mm (length) × 6.86 mm (width) × 2.54 mm (thickness). The process parameters reported are 750 W laser power, 12.7 µm raster scan spacing, and 600 mm/s raster scan speed. The raster scanning speed (v R ) is related to the linear velocity (v) of the laser in the X direction by v = scan spacing×v R 2×width . Figure 5. Longitudinal cross section of a CMSX-4 specimen fabricated with 750 W laser power and scan speed of 600 mm/s showing increase in melt pool depth along the length of the specimen (Reproduced with permission from the authors of [51].). Figure 6 shows the comparison between the analytically calculated and experimentally observed melt pool depth for the CMSX-4 specimen as illustrated in Figure 5. The maximum absolute relative error between the experimental and analytical melt pool depths is~8% at the start of the scan period, with the mean absolute error during the 30 s period being 2.71%. The results also show the transient nature of the melt pool, with the depth continuously increasing as a function of time with a fixed set of process parameters.

Hybrid Model
The validation results show that Eagar-Tsai's model can provide excellent low-cost estimates of melt pool dimensions for single-track and single-layer AM depositions. The goal of this work is to obtain model-based estimates of the required temporal variations in process parameters for achieving target melt pool dimensions during the build process. In order to achieve this, an optimization problem that finds process parameters as a function of time is solved with an objective of minimizing the deviation from the desired melt pool dimensions. As the true functional form that relates the process parameters to the melt pool dimensions is unknown, the information regarding gradients is not readily accessible, therefore the optimization problem is in a black-box setting. Several such black-box optimization techniques are well-studied by different researchers [52], e.g., stochastic process-based approaches (kriging methods) [53], evolutionary algorithms [54], trust-region based algorithms [55], and random search [56]. These approaches typically involve iterative sampling of the objective function in the search space. Such sampling often proves to be prohibitively expensive, particularly in high dimensions, when the involved process model is computationally very expensive.
To address this challenge, in the present work, an application of gradient-free optimization under budget limitations is implemented by solving the melt pool geometry control problem using BO that employs computationally inexpensive GP surrogate models learnt via sparse sampling of the space of process parameters (e.g., P and v). The details of the mathematical formulations of GP and BO can be found in Appendix A. Although all optimization demonstrations in this paper involve the Eagar-Tsai's model as the process model, it can very well be replaced by any other process models of users' choice, without any fundamental change in the optimization algorithm.

Results and Discussion
This section presents the results and discusses the significant findings therefrom.

Prediction of Melt Pool Dimensions Using GP
The melt pool dimensions in the L-PBF AM process continuously increase due to progressive heating of the specimen caused by the layer-wise fabrication as evinced by simulation results [57] as well as experimental observation [47]. Figure 7 shows the variation in melt pool dimensions as a function of time for a CMSX-4 specimen fabricated using a laser power of 900 W and linear scan velocity of 0.25 mm/s. The melt pool depth increases from 1.814 mm at t = 1 s to 3 mm at t = 20 s, while the melt pool width increases from 3.82 mm to 6.11 mm. The percentage increase in melt pool depth is~65%, and the increase in width is~59%. The melt pool reaches a steady state at~20 s. If the melt pool dimensions need to be maintained at desired values, the process parameters should be adjusted during this transience period. However, performing simulations for a range of process parameters over the entire deposition process is computationally expensive when coupled with an optimization framework. In applications where data is scarce, particularly when each simulation can be potentially computationally expensive, it is often infeasible to have a large number of training points. Moreover, inherent uncertainties in the physical process and the simulation models (e.g., uncertainties in the thermophysical properties of the material) mandate the need of uncertainty quantification with the predictions [33,58]. Surrogates, like GPs, can be very useful in such cases in order to have probabilistic estimates of the quantity of interest with limited datasets, whereby the posterior mean and variance can not only be informative for prediction at unknown input locations, but the information can also be used for budget-constrained optimization, as described in the following subsection. A surrogate modeling and optimization technique as discussed in Appendix A is adopted in this work to find the model-based estimates of the process parameters required for controlling the melt pool dimensions.
In order to build the surrogate using GPs, two-hundred (200) Latin hypercube sampling (LHS) simulations are performed to obtain the values of melt pool depth and width at each time instant for combinations of P and v in the range of 300 to 1200 W and 0.5 to 2.5 mm/s, respectively. This parameter range is chosen to avoid any keyhole mode of melt pool formation in CMSX-4 where the Eagar-Tsai's model shows limited accuracy [59,60]. The cross section of melt pools created in conduction mode is generally semicircular, as predicted by Eagar and Tsai's conduction mode model [35]. LHS is selected as this is one of the most commonly used statistical methods for DoE. It allows for a good spread of the initial DoE over the design region with limited iterations due to its high sampling efficiency [61]. For every time instant under consideration, each DoE point comprises a combination of process parameters (P, v) and its corresponding melt pool depth and width values.
The effect of training data size on the regression performance of the GP surrogate in predicting melt pool depth and width is of critical interest. Out of the 200 initial LHS points, a test set of randomly selected 100 points is set aside for testing the regression performance of the GP models. Different GPs are trained at different time instants with randomly selected samples from the training set. The training data size is varied from 10 to 100, in steps of 10 samples. Ten random selections of the training data is chosen for each training size, in order to find the average behavior of the regression performance for each training size. The prediction performances are tested on the set aside 100 samples for each trained GP model.
The metric for gauging the prediction performance is chosen as the Relative Squared Error (RSE) ||ŷ−y * || 2 ||y * || 2 , whereŷ is the prediction and y * is the true value (i.e., ground truth as obtained by the Eagar-Tsai's model). A lower RSE suggests that the most likely prediction (mean of the predicted posterior Gaussian distribution) of the depth/width in the test set matches closely with the true value. Figure 8 shows that the RSE is quite low for all the time steps, while there is a slight increase in RSE from 2 s to 20 s for both depth and width. RSEs for the melt pool depth and width prediction show a sharp drop from training size of 10 to 20 samples, after which the RSE saturates for almost all the higher training sizes. This indicates that a training size of at least 20 LHS DoE points is in general sufficient for a satisfactory prediction of the melt pool depth and width in the selected process parameter space. This provides an estimate of the number initial DoE points required for learning the surrogate model for the subsequent BO steps at each time instant.   Figure 10 for the melt pool depth and width predictions at 20 s, with a training set size of 20 samples. The R 2 score decreases to 0.82 for depth prediction and 0.88 for width prediction, which is also explained by the higher RSE in prediction at 20 s (Figure 8). The variability in the predicted values from the true values is explained by the wider ±2σ band at the test locations. The lower R 2 score at 20 s as compared to 2 s can be attributed to the higher variability in the steady state values of depth and width as compared to the initial stages as a function of the process parameters. From the perspective of applicability of the surrogate predictions, GPs provide the end user with not only a computationally inexpensive way of predicting the melt pool dimensions at different operating conditions for which experiments and simulations are not performed, but also with an estimate of uncertainty quantification (UQ) for predictions from the variance associated at the query points.

Optimization of Process Parameters for Melt Pool Dimension Control: Objective Function
As discussed in Sections 1 and 2, appropriate control of the process parameters, e.g., P and v are required for maintaining desired melt pool dimensions during the deposition process. The goal of controlling the melt pool dimensions is formulated as an objective function that needs to be maximized in the setting of BO, as described in Appendix A.3. In order to pose the melt pool control problem in an optimization setting, an appropriately defined objective function (J) is to be maximized at every discrete time instant of control, where d and W denote the depth and width, respectively, at a particular time instant, whereas d * and w * represent the desired depth and width, respectively, during the deposition process. P indicates the power at time instants of control, whereas P max and P min denote the maximum and minimum values of the range of laser power in the space of process parameters. Therefore, |P−P min | |P max −P min | denotes the normalized power input at a particular time instant, and incorporating it into the objective functional serves as a penalty term for high power values. This is introduced since it is desired to achieve the controlled process along with avoiding processing conditions involving very high levels of power. c 1 , c 2 , and c 3 in Equation (3) denote the relative weights of the components of the objective function that controls the depth, width, and power, respectively. By varying the relative values of c 1 , c 2 , and c 3 , it is possible to preferentially weigh the objective function to meet the requirements. Formulating the objective functional J is a key step in the optimization process, and the optimization routine should be designed in such a way that undesirable process parameters result in low values of objectives. This formulation conforms to the mathematical characteristics of the Matérn covariance function [15] used in this work, which assumes local smoothness of the inputs, so that input process parameters that are close together in the (P, v) space are expected to have similar objective values. As undesirable parameter combinations have low objectives, points in the close vicinity of them will have low acquisition potential (see Appendix A) during the optimization steps, and therefore will have lower priority for the incumbent selection.

Optimization Routine
A sequential global-local optimization thread [62] is employed in this work for the selection of process parameters for controlling the melt pool geometry. The terms "global" and "local" correspond to the search spaces with respect to which BO is performed. In the global optimization thread, a potential optimal process parameter combination is selected, around which a refinement is made during the local optimization thread for the final selection of the optimal process parameter combination. The outline of the process is described in the flowchart as depicted in Figure 11. Details of the optimization algorithm can be found in Appendix B. Entities capped with "tilde" (∼) correspond to parameters of the local optimization thread. The control problem is designed to maintain the melt pool depth (d * ) at 2.5 mm and width (w * ) at 5.0 mm. The initial global GPs are learnt with 20 N init LHS points in the (P, v) space, with P ∈ [300, 1200] W and v ∈ [0.5, 2.5] mm/s. The ranges of the process parameters and the controlled geometrical parameters chosen in this problem are often dictated by the experimental requirements and constraints, and can be modified as needed. The optimization problem over the entire duration of the L-PBF process can be considered as a sequence of optimization problems performed at some discrete intervals of time, which motivates the choice of learning separate surrogates for each discrete time instant for which the melt pool needs to be controlled.
With the initial trained GP from 20 LHS samples for each time instant, low-cost surrogate predictions are made over the large global search space X star , consisting of 2000 uniformly chosen points. This methodology forms the key to the surrogate-based optimization processes: predictions are made over an extensive search space by employing a low-cost surrogate by avoiding functional evaluations (physics-based simulations in this case) at the search points. Based on the mean-variance trade-off of the GP's posterior predictive distribution, the acquisition function EI (See Appendix A) guides the search iteratively to subsequent optimization points until the computational budget is depleted. N Global Iter = 10 optimization steps are chosen for the global thread, after which the optimal process parameter combination with the highest objective value is selected.
For all the time instants of control, it was possible to obtain optimal process parameters within the global thread that can yield the melt pool depth and width close to the desired values.
However, for finer adjustments in P and v values in order to achieve the micron-level control in the melt pool geometry, the local optimization thread is executed. Here, the local search spaceX star is built aroundx * , the optimum with the highest objective value from the global thread, in the following way; P is varied within ±20 W and v within ±0.15 mm/s of the corresponding values inx * . 20Ñ init LHS points are chosen within the limits ofX star , the optimization routine is performed within this local search space forÑ Local Iter = 10 steps. This yields the refinement of the process parameter values that result in melt pool dimensions very close to the ones desired. The number of optimization iterations in both the local and global threads (i.e., N Global Iter andÑ Local Iter ) are kept as design parameters in this paper, which are expected to be driven by computational budget for practical applications. The scaling factors for the objective function components are chosen as c 1 = 1, c 2 = 0.1, and c 3 = 0.1, which are based on the greater relative importance of maintaining the depth as close to d * as possible during the build process, so that a uniform deposit is maintained. Figure 12 shows the performance of the surrogate based optimization method. Figure 12a,b shows the variation of the controlled melt pool depth and width from the desired d * and w * values as a function of time. It is seen that the maximum variation in the meltpool depth is 1.1 µm which is 0.04% of the desired depth, while that for width is 173 µm, which is 3.46% of the desired width. The process parameters change from P = 948.07 W, v = 0.40 mm/s at t = 2 s to P = 737.81 W, v = 0.37 mm/s at t = 20 s (Figure 12c,d). It is to be noted that penalization of P in the objective function allows us to have a smoothly variation of P over 20 s with lower P values, which changes by~22% of the starting value at t = 2 s. The maximum change in v is~27%. It is possible to have process parameters with higher P values (and higher v values) that control the melt pool, but those conditions are avoided by P penalization in the objective function. Similar formulation of the objective function can be pursued by v penalization, if smother v transition is desired.

Validation with Experimental Results of Melt Pool Depth Control
The model-based control strategy is validated with experimental results reported in the open literature [63]. Figure 13a shows the longitudinal cross section of an L-PBF processed René 80 specimen fabricated using a powder layer thickness of 1.4 mm on a substrate of dimensions 35.56 mm × 6.86 mm × 2.54 mm. The experiment was carried out with the raster scanning speed of 450 mm/s and scan spacing of 25.4 µm. In reality, experiments are extremely challenging to perform with melt pool depth as a controlled variable as there exists no easy way of measuring the melt pool depth in situ. During the conduction mode, the surface temperature control of melt pools was found to correlate well with the melt pool depth experimentally by Bansal et al. [63]. The observation was also made by Raghavan et al. in E-PBF of Inconel 718 [64]. The mean value of the melt pool depth during the control period is around 1200 µm from experiments. Accordingly, d * is set at 1200 µm for the surrogate-based optimization routine involving only P as outlined in the experiments.
The melt pool simulation is performed by employing the Eagar-Tsai's model, with the thermo-physical properties of solid René 80: k = 24.56 W/(m·K), ρ = 7604 kg/m 3 , c = 600 J/(kg·K), and the liquidus temperature T L = 1607 K obtained using the software JMatPro [65], with the alloy composition of René 80 provided by the vacuum alloy product catalog of Cannon Muskegon [66]. The objective function involves two components in this case: with c 1 = 1 and c 2 = 0.1. Figure 13 shows that the predicted melt pool depth with the surrogate assisted control scheme remains very close to the desired d * . A variation of ∼200 µm is observed in the experimental results, whereby the melt pool depth shows a slightly increasing trend towards the end of the control process, reflected by the slight increase in power input from 10 s till 14 s. The trend of the surrogate predicted controlled power input matches closely with the experimental results till 8 s, which is explained by the fact that for a constant laser speed, the power input should decrease as a function of time in order to maintain a constant melt pool depth as the specimen gains thermal energy continuously. Nonetheless, the validation results indicate the efficacy of the surrogate-based optimization routine in predicting process parameters as a function of time for achieving melt pool dimension control during the deposition process.

Summary and Conclusions
This paper develops a novel hybrid methodology for control of melt pool geometry in AM in the framework of ML-assisted modeling and optimization. The continuous changes in the melt pool geometry are predicted by a low-cost GP surrogate developed using an experimentally validated analytical 3D model. The uncertainties in the predictions for the melt pool geometry are quantified using GPs. Reliable estimates of the optimal process parameter evolution are obtained through active learning via BO by devising appropriate objective functions that quantify the control requirements. The methodology provides an estimation of process parameter variation during the AM deposition process in order to maintain a desired melt pool geometry under continuously varying thermal conditions.
Being data-driven, the reliability of the optimized process parameters obtained from the algorithm is based on the accuracy of the underlying physical model in predicting the melt pool geometry in the range of process parameters considered. However, with a high-fidelity model, the computational cost involved in the optimization process can be significantly high. For example, the cost of running a Netfabb [67] DED model for a single-track and single-layer process on a CMSX-4 specimen as described in this paper is ∼10 times as expensive as Eagar-Tsai's model at t = 2 s, and ∼150 times at t = 10 s. Although such high-cost simulation is prohibitive in practical applications, in the future, the methodology will be enhanced by using information from different fidelities [68] of AM computational models to develop a multifidelity modeling and optimization framework for prediction and control of melt pool geometries. For example, a Netfabb model can be combined with inexpensive observations from Eagar-Tsai's model to develop a two-fidelity GP that can be used for melt pool geometry prediction and optimization with an expected reduction in computational cost. Additional investigations are also planned as summarized below.

1.
Usage of process parameters (other than P and v) such as scan spacing (i.e., hatch spacing), scan pattern (i.e., hatch pattern), build plate temperature, and powder layer thickness (for powder bed AM) and powder feed rate (for directed energy AM) as control inputs.

2.
Incorporation of design constraints in the surrogate assisted modeling framework, e.g., tackling harder problems whereby a melt pool geometry needs to be controlled, along with the final microstructure such as columnar grains in CMSX-4 .

3.
Development of heterogeneous design spaces in formulating multifidelity modeling framework, which can be useful in catering to optimization problems where the varied levels of fidelities have different input spaces. As an example, a HF L-PBF process model can have several process parameters, such as powder distribution properties, hatch spacing, scan strategy, etc., apart from P and v of the laser (which are the only process parameters for a LF model like the Eagr-Tsai's model), which can potentially affect the build characteristics. Optimization in such a framework can be possible via heterogeneous transfer learning [69] to learn from a common subspace of the inputs.

4.
Integration of the surrogate modeling framework developed in this work with microstructure prediction framework using an open source SPPARKS code [70] to optimize the process parameters for obtaining a desired microstructure.

5.
Implementation of the surrogate modeling framework for checking its efficacy towards mitigating unintended melt pool behavior such as keyholing effect. However, in such a case, the computational model needs to be a high fidelity one capable of predicting such a behavior.
It is envisioned that the newly developed framework can be implemented in developing feedback strategies for melt pool control during AM processes with shorter response time due to improved prediction capability [71]. Remark A1. The main advantage in having a probabilistic prediction method is that it naturally provides a measure of uncertainty quantification (UQ) associated with these predictions through the variance of the distribution. Moreover, having a probabilistic estimate instead of a fixed estimate provides the end user with a confidence level of the predictions.
For a finite subcollection {x 1 , x 2 , · · · , x N } of the random input x, the corresponding objective values are assumed to have a multivariate jointly Gaussian distribution: , and E(·) indicates the expectation of a random variable.
Let D trn = {(x trn i , y trn i )}, i = 1, · · · , N, be the training set that are available for model formulation. For a noisy regression model, it is assumed that where the additive noise ε i are independent and identically distributed (iid) zero-mean Gaussian random variables, i.e., ε i ∼ N (0, σ 2 ). By incorporating the noise term, the joint distribution of the observed values and the functional values at the test locations are where D tst = {(x tst , y tst )} is the test set in which y tst is online observed data and x tst is the unknown variables to be estimated. Therefore, by the property of multivariate normal distributions, the conditional distribution of the function values at the test location is Gaussian [72]. In particular, Thus, the algorithm predicts the mean and covariance for the posterior distribution that models the output at every test data point.
are aimed at finding a trade-off between exploration and exploitation in possibly noisy settings [52,75,76], which facilitates a balance between global search and local optimization through acquisition functions. Usage of GPs in the surrogate modeling framework for BO often leads to closed-form analytical expressions for acquisition functions, which are inexpensive to compute.
In the above formulation, the acquisition functions are designed in the following way. The potential of performance improvement is driven by the predicted mean function (which is in the category of exploitation), whereas the uncertainty prediction is manifested by regions of high variance (which is in the category of exploration). A trade-off between these two requirements is achieved by the acquisition functions iteratively, as a sequential optimization process. One commonly used acquisition function in Bayesian optimization is Expected Improvement (EI), which is employed in this work. According to the formulation of Mockus et al. [77] and Jones et al. [52], the EI acquisition function can be written as and x + = argmax x i ∈x 1:k f (x i ) is the input corresponding to the maximum functional value sampled until iteration k; the parameter ξ > 0 controls the trade-off between exploration and exploitation [31]; µ(x) and σ(x) are the mean and variance, respectively, predicted by the GPs for the an input point x; and φ(·) and Φ(·) are the probability distribution function (PDF) and cumulative distribution function (CDF) of the standard normal distribution, respectively. If the objective function is noise corrupted, instead of using the best observation, the point with the highest expected value is defined as the incumbent, i.e., f (x + ) replaced by µ + , which is defined as The optimization algorithm proceeds sequentially by samplingx = argmax x EI(x) at every step of the iteration process to add on to the dataset, after which the GP model is retrained with the new data set to predict the acquisition potential for the next iterative step. This process continues until an optimum is reached, or the computational budget is extinguished. As the acquisition potential is predicted over the entire search space by the surrogate, BO can achieve fast predictions without a lot of function calls in the search space (i.e., without having to run the simulations to obtain the objectives at all the search locations), which, otherwise, might be computationally infeasible when the search space is high-dimensional and the simulations are expensive. Figure A1 demonstrates the optimization process through a simple problem of maximizing a black-box function f (x) in the input space of [−2, 2] (denoted by solid blue line). The process is started with an initial design of experiments (DoE) of 4 points (filled blue circles) chosen in the search domain by Latin hypercube sampling (LHS) [61], which allows the initial selection of points to be well-spread in the search space. Figure A1a shows the GP prediction in the initial stage, which is characterized by the mean function (red dotted line), and the variance associated with the prediction (denoted by the orange band). The variance is high globally at the initial stage, because only four points are sampled on the true function. The maximizer of the acquisition function EI(x) shows the location to sample for the next point at every iteration (filled red circles). Figure A1b shows that the predicted mean of the function in Figure A1a drives the selection of the first optimization point. It is interesting to note that in Figure A1c, although the true maximizer of f (x) is already sampled, EI(x) points towards a region where the variance is high. This is because of the exploration-exploitation trade-off which allows the algorithm to explore the search space instead of greedily searching for the optimum, a property that helps the algorithm avoid being stuck at a local optima. Eventually after seven optimization iterations, the routine converges at the true optima, as shown in Figure A1d. Moreover, the mean function predicted at the last stage is very close to the true function, and the uncertainty band also reduces globally, which indicates that the algorithm not only finds the maximizer of the function, but in doing so it also learns a pretty accurate surrogate model of the function, which can now be used as a low-cost approximation of the true function.

Appendix B. Optimization Algorithm
Algorithm A1 shows the surrogate-based optimization routine that is employed in this paper for solving the melt pool geometry control problem, which is elucidated in Section 3.3. The global thread picks up an optimal pointx * , which is refined in the local thread to give x * , as the final optimal point chosen by the algorithm.

Algorithm A1 Optimization Routine Global Thread
Require: Start with surrogate GP for the objective functional J with N init LHS initialization for k = 1 to N Global Iter optimization steps do -Use the trained GP to predict the posterior distribution of the objective function in the global search space X star -Sample an incumbent optimized process parameter from X star as x * k = argmax x EI(x) -Compute the objective function J k at the chosen x * k -Add the pair (x * k , J k ) to the surrogate model's input and output sets respectively -Retrain the surrogate model end for -Selectx * = argmax x k J -SelectX star appropriately aroundx *

Local Thread
Require: Start with surrogate local GP for the objective functional J withÑ init LHS initializations for k = 1 toÑ Local Iter optimization steps do -Use the trained local GP to predict the posterior distribution of the objective function in the local search spaceX star -Sample the optimized process parameter fromX star asx * k = argmaxx EI(x) -Compute the objective function J k at the chosenx * k -Add the pair (x * k , J k ) to the local surrogate model's input and output sets respectively -Retrain the local surrogate model end for -Select x * = argmaxx k J