Abstract
In selective laser melting (SLM), the relationships between process parameters and surface quality are highly nonlinear, while different quality indicators are often mutually conflicting, making multi-objective optimization under limited experimental budgets challenging. To address this issue, this study proposes a direct Gaussian process (GP)-based multi-objective Bayesian optimization framework for the efficient optimization of SLM process parameters. Targeted improvements in surrogate modeling and acquisition function design enable effective characterization of coupled multi-objective responses with a limited number of evaluations. Surface roughness (Sa), defect area ratio, and macroscopic warpage root mean square were selected as optimization objectives based on laser confocal microscopy measurements, while laser power, scanning speed, and hatch spacing were treated as design variables. Comparative studies with the conventional GP-expected hypervolume improvement method and the evolutionary algorithm Non-dominated Sorting Genetic Algorithm II (NSGA-II) demonstrated that the proposed approach achieved faster convergence and superior Pareto front approximation under the same evaluation budget. From an engineering perspective, the proposed framework reduces reliance on extensive experiments and enables the identification of practically feasible parameter combinations, providing an efficient and transferable solution for intelligent SLM process optimization.
Introduction
Selective laser melting (SLM), as one of the most representative technologies in metal additive manufacturing, has attracted extensive attention owing to its capability for fabricating complex geometries, achieving high dimensional accuracy, and enabling near-net-shape production.1–3 Compared with conventional subtractive manufacturing, SLM realizes the integrated formation of complex components through a layer-by-layer melting and deposition mechanism, demonstrating significant potential in aerospace, biomedical, and high-performance engineering applications.4–6 However, the surface quality and process stability of SLM-fabricated parts are highly sensitive to processing parameters; slight variations in laser power, scanning speed, layer thickness, and hatch spacing can markedly influence the final build quality.7,8
Processing parameters regulate melt pool dynamics, energy input, and solidification behavior, thereby affecting critical quality metrics such as surface roughness, porosity formation, residual stress accumulation, and macroscopic warpage deformation.9,10 Notably, surface roughness, defect level, and warpage often exhibit trade-off relationships under different parameter combinations. 11 For example, increasing energy density may reduce lack-of-fusion defects while simultaneously intensifying balling phenomena or deteriorating surface roughness. 12 Consequently, achieving coordinated optimization of SLM process parameters under multiple, potentially conflicting quality objectives remains a significant challenge.
Traditional parameter optimization approaches primarily rely on design of experiments, response surface methodology (RSM), or empirical trial-and-error strategies.13,14 Although these methods provide a statistical basis for analysis, they generally assume relatively simple response relationships and require substantial experimental data, limiting their applicability in high-dimensional and strongly nonlinear parameter spaces. 15
In recent years, neural networks, particle swarm optimization, and multi-objective evolutionary algorithms (e.g., NSGA-II and its variants) have been widely employed for Pareto optimization of SLM process parameters. Zhang et al. 16 developed a backpropagation neural network model with laser power, scanning speed, and hatch spacing as inputs and surface roughness as the output, providing a predictive basis for parameter optimization. Theeda et al. 17 constructed a unified artificial neural network to simultaneously predict multiple performance indicators, including surface roughness, relative density, microhardness, and dimensional accuracy. Shirmohammadi et al. 18 established a multilayer perceptron model combined with particle swarm optimization for parameter tuning while introducing RSM-based quadratic regression to analyze parameter effects and determine optimal combinations for surface roughness control. Murat et al. 19 built an RSM-based surrogate model linking process parameters to tensile properties and integrated it with particle swarm optimization to reduce experimental effort while achieving effective parameter prediction. Detwiler et al. 20 developed decision tree and neural network models for the high-accuracy prediction of surface morphology and revealed the critical roles of thermocapillary convection and melt track overlap in surface topology formation. Cao et al. 21 constructed a literature-based database and applied a clustered ensemble regression model to enhance prediction accuracy, combined with NSGA-II for strength–ductility multi-objective optimization, and employed SHapley Additive exPlanations analysis to interpret the synergistic regulation mechanisms of process parameters on microstructure and properties in L-PBF Ti6Al4V. Li et al. 22 developed ensemble surrogate models to characterize the relationships among process parameters, energy consumption, tensile strength, and surface roughness and conducted NSGA-II-based multi-objective optimization with experimental validation. Gao et al. 23 proposed a Gated Recurrent Unit neural network model based on interlayer infrared radiation intensity sequences to predict thermal accumulation behavior in Laser Powder Bed Fusion (LPBF) and improved density, surface quality, and residual stress through parameter optimization.
However, population-based algorithms typically require a large number of function evaluations to achieve convergence. When each evaluation corresponds to a physical experiment or high-fidelity numerical simulation, the associated computational and experimental costs become prohibitive for practical engineering applications.24,25
Bayesian optimization (BO), as an efficient strategy for optimizing expensive black-box functions, has received increasing attention in recent years. 26 By constructing a Gaussian process (GP) probabilistic surrogate model, 27 BO captures complex nonlinear mappings and guides sequential sampling through uncertainty quantification. Multi-objective Bayesian optimization (MOBO) further extends this framework by incorporating acquisition strategies such as expected hypervolume improvement (EHVI) to efficiently address conflicting objectives.28,29
In the field of additive manufacturing, several studies have explored the application of BO and surrogate-based methods for process parameter optimization. Gnanasambandam et al. 30 proposed a deep GP-based BO approach combined with stochastic imputation inference and bootstrap aggregation to model and optimize nonstationary process responses. Heddar et al. 31 developed a scalable multi-objective optimization framework using space-filling designs to construct high-fidelity GP surrogate models, integrating BO with evolutionary algorithms to jointly optimize porosity and surface roughness. Chepiga et al. 32 employed GP-based surrogate modeling with MOBO to achieve Pareto optimization of hardness and porosity in laser powder bed fusion under limited experimental data. Ding et al. 33 proposed an experimental data-driven MOBO framework utilizing EHVI to coordinate printing quality and efficiency while reducing experimental costs.
Nevertheless, in SLM scenarios, existing MOBO methods are commonly built upon standard EHVI frameworks or indirect surrogate modeling strategies, which may suffer from limited sampling efficiency and insufficient capability to capture strongly coupled multi-objective responses under restricted evaluation budgets.34,35
To address these limitations, this study proposes a direct GP-based MOBO framework for efficient SLM process parameter optimization. Compared with conventional MOBO implementations, targeted improvements are introduced in surrogate model construction and acquisition function design to enhance sampling diversity and convergence efficiency. Based on laser confocal microscopy measurements, surface roughness (Sa), defect area ratio, and macroscopic warpage root mean square (RMS) are selected as optimization objectives, while laser power, scanning speed, and hatch spacing are treated as design variables to establish an engineering-representative optimization space.
Comparative studies with the traditional GP-EHVI method and the multi-objective evolutionary algorithm NSGA-II demonstrate that, under the same evaluation budget, the proposed method achieves faster convergence and superior approximation of the Pareto front. From an engineering perspective, the framework significantly reduces reliance on large-scale experiments or computationally expensive simulations, facilitating the identification of practically feasible parameter combinations under limited experimental conditions. This work provides an efficient and generalizable pathway for the intelligent optimization of SLM process parameters and offers new insights into data-driven optimization methodologies in metal additive manufacturing.
Material and Methods
Experimental data collection
Surface measurement
This study focuses on the microscale surface morphology analysis of metal components fabricated by SLM. Owing to the structural complexity and abundance of microfeatures inherent to SLM-processed surfaces, high-precision measurement techniques are required. Therefore, a laser confocal microscope was employed for surface characterization in this work.
Laser confocal microscopy offers the advantages of high spatial resolution and noncontact measurement, enabling accurate capture of fine surface features such as surface roughness, layer-wise textures, and porosity. Its efficient three-dimensional imaging capability and superior depth resolution make it particularly suitable for observing the complex layered structures and microscale features generated during the SLM process. In addition, laser confocal microscopy provides images with high contrast and clarity and supports automated analysis, thereby significantly enhancing both measurement accuracy and efficiency.
The measured surface images are shown in Figure 1a. As can be observed, pronounced vertical textures are present on the surface, indicating relatively poor surface quality accompanied by numerous defects, such as spherical spatters. Figure 1b presents the height map of the sample surface, where a large number of microscale protrusions can be clearly identified. These protrusions are likely attributable to spherical spatter generated during the printing process as a result of excessive laser power. Moreover, the surface exhibits highly pronounced texture characteristics.

Real surface image and three-dimensional surface height map.
Data acquisition
The surface quality of components fabricated by SLM is governed by the coupled interactions among the laser, powder bed, and molten pool, resulting in surface characteristics that vary across multiple spatial scales. To ensure a consistent and reliable data basis for subsequent MOBO, this study focuses on three representative surface quality indicators that capture the dominant surface characteristics relevant to process optimization.
At the microscale and local defect level, the arithmetic mean surface roughness Sa and the defect area ratio are adopted as the primary evaluation metrics. The roughness parameter Sa provides an overall description of surface height variations, while the defect area ratio quantitatively characterizes the extent of surface defects such as spatter particles and balling. Both metrics exhibit clear physical meaning and favorable stability under different process parameter conditions, making them suitable as objective responses for multi-objective modeling.
At the macroscopic scale, global warpage deformation induced by thermal stress accumulation during the SLM process was considered. A low-frequency component extraction method was employed to isolate macroscopic warpage from the complex surface morphology. The RMS value of the extracted warpage surface was then used as the quantitative indicator of global deformation, which is defined as follows:
Based on the above characterization results of surface morphology and macroscopic deformation, Sa, defect area ratio, and warpage RMS were selected as the optimization objectives in the MOBO framework. Table 1 summarizes the experimental dataset used for GP-based MOBO. Unlike the previous simplified description involving only five representative samples, the actual dataset consists of 20 distinct process parameter combinations, with three repeated samples collected for each combination, resulting in a total of 60 experimental samples.
Data Acquisition Results of the Five Samples
HS, hatch spacing; LP, laser power; RMS, root mean square; SS, scanning speed.
The dataset includes the measured surface quality indicators, namely arithmetic mean roughness (Sa), defect area percentage, and macroscopic warpage RMS, together with the corresponding process parameters, including laser power (LP), scanning speed (SS), and hatch spacing (HS). The repeated measurements under each parameter combination were introduced to reduce experimental uncertainty and improve the statistical robustness of surrogate model training. The resulting 60-sample dataset spans multiple regions of the three-dimensional process parameter space and provides sufficient information for GP hyperparameter estimation and posterior uncertainty modeling.
To further evaluate surrogate model reliability, five-fold cross-validation was performed during GP training, and the corresponding prediction errors were monitored throughout the optimization process. Therefore, the collected dataset provides the quantitative basis for subsequent MOBO and enables analysis of the nonlinear relationships between process parameters and surface quality characteristics.
Algorithms and models
Optimization problem formulation and objective function definition
In this study, 316 stainless steel powder was selected as the build material, and an optimization investigation was conducted focusing on key process parameters that exhibit high sensitivity to energy input and forming quality during the SLM process. Considering the processing capability of the employed equipment, the stable forming window, prior experimental experience, and relevant literature, LP, SS, and HS were identified as the primary optimization variables. The corresponding parameter ranges are summarized in Table 2.
Process Parameters and Their Ranges
In this work, the layer thickness (LT) was fixed at 30 µm to ensure process stability and adequate interlayer metallurgical bonding. On this basis, the volumetric energy density (VED) was introduced as a comprehensive indicator to characterize the laser energy absorbed per unit volume of material. The VED is defined as follows:
Under the specified ranges of process parameters and the fixed LT condition, the resulting VED spans a typical process regime from low to high energy input. This range is capable of comprehensively capturing variations in melt pool behavior, forming quality, and defect evolution under different energy density levels, thereby providing a sufficiently representative and well-defined design space for establishing nonlinear mappings between process parameters and performance responses in subsequent multi-objective optimization.
To address the requirement for the simultaneous optimization of multiple quality-related metrics in the SLM process, the process parameter optimization task was formulated as a surrogate model-based multi-objective minimization problem. Within this modeling framework, the optimization variables, objective functions, and constraints are consistently defined according to the actual process parameter space and experimental data, and an efficient solution was achieved through a MOBO strategy.
The process parameter vector is defined as follows:
With respect to objective function construction, three key indicators that comprehensively reflect forming quality and structural stability were selected as optimization objectives: surface roughness Sa, defect area fraction (defect percentage), and warp RMS. All these metrics belong to the category of “smaller-is-better” performance indicators and can therefore be uniformly formulated as a multi-objective minimization problem:
It is important to note that these quality metrics—surface morphology features, thermal deformation indicators, and defect characterization parameters—cannot typically be described by explicit analytical models, as they exhibit strong nonlinearity, high coupling, and noise sensitivity with respect to the process parameters. Consequently, a data-driven black-box modeling strategy was adopted, in which each objective is approximated by a GP regression model that provides predictive means and uncertainty estimates under small-sample conditions. This black-box formulation naturally aligns with the BO framework described below.
Multi-objective Bayesian optimization
MOBO is a class of global optimization methods designed for high-dimensional problems involving expensive evaluations and conflicting objectives. It has been widely applied in experiment-driven engineering optimization and manufacturing process parameter design. The core idea of MOBO is to construct probabilistic surrogate models of objective functions under a limited experimental budget, thereby balancing exploration and exploitation to progressively approximate the true Pareto-optimal frontier.
A general multi-objective optimization problem can be formulated as follows:
BO addresses this challenge by modeling each objective function using a GP surrogate, which provides a probabilistic description of the objective function uncertainty. For the m-th objective, the GP prior can be expressed as follows:
In the multi-objective setting, MOBO typically adopts the hypervolume (HV) indicator as a unified optimization criterion. Given a reference point
Based on the HV indicator, the EHVI is introduced as the acquisition function. EHVI is defined as the expected increase in HV contributed by adding a new sample point x under the GP posterior distribution:
From a probabilistic perspective, EHVI encourages sampling in regions with high potential near the current Pareto front (exploitation) while also promoting exploration in regions with large predictive uncertainty, thereby enabling adaptive guidance of the multi-objective optimization process.
This uncertainty information plays a critical role in the construction of acquisition functions. In conventional pipelines, variance information is either discarded after EHVI integration or buried inside the Monte Carlo (MC) noise; here it is retained as an accessible quantity that can be directly fed into the acquisition function and the batch selector, forming the basis for the direct coupling channels described above.
Based on the aforementioned theoretical foundation, a direct GP–EHVI coupling framework was developed in this study (as illustrated in Fig. 2). Unlike conventional MOBO, where the GP merely feeds the predictive distribution into the EHVI estimator as a passive input, the proposed framework establishes three explicit coupling channels: (C1) the GP predictive standard deviation is analytically added to the EHVI scalar as an explicit exploration bonus; (C2) the GP predictive mean directly drives objective-space diversity assessment during batch selection; and (C3) the trade-off between exploitation and exploration is adaptively regulated according to the HV convergence rate. These channels are detailed in “Construction of multi-objective acquisition functions and optimization strategy” section.

Direct Gaussian process-based multi-objective Bayesian optimization framework.
The optimization starts from a limited set of initial samples. In each iteration, the surrogate models are updated, the acquisition function is optimized to select new candidates, and the Pareto front is accordingly refined. The process continues until a predefined evaluation budget or convergence condition is met. Under limited evaluations, the proposed method demonstrates improved convergence efficiency and well-distributed Pareto solutions in comparative studies, highlighting its effectiveness for multi-objective parameter optimization in complex manufacturing processes.
Construction of multi-objective acquisition functions and optimization strategy
To simultaneously balance exploration and exploitation, the EHVI is adopted as the acquisition function. Because an analytical expression for EHVI is generally unavailable in the multi-objective case, a MC-based numerical approximation is employed. The reference point is dynamically updated based on the worst observed performance in the current sample set, ensuring that the HV calculation always reflects the true potential for improvement.
However, in standard GP-EHVI, the predictive uncertainty
To further enhance efficiency under limited evaluations, a batch sampling strategy with diversity enhancement is incorporated. Conventional batch strategies typically impose diversity constraints solely in the decision space or require true objective values to compute objective-space distances—information unavailable before evaluation. To overcome this limitation, the second coupling channel (C2) utilizes the GP predictive mean to guide objective-space diversity. The first point in each batch is selected by maximizing the uncertainty-augmented EHVI:
For subsequent points (
The use of
Finally, the diversity weight
To clarify the structural distinction between the proposed framework and conventional MOBO, Table 3 provides a step-by-step algorithmic comparison. While standard GP-EHVI maintains a unidirectional information flow (GP → EHVI) with fixed or preset diversity weights, the direct GP-EHVI framework establishes bidirectional coupling through the three channels defined above: C1 injects GP predictive uncertainty directly into the acquisition scalar, C2 employs GP predictive means for objective-space diversity assessment without requiring true evaluations, and C3 enables adaptive regulation of coupling strength based on optimization-history feedback. These differences are instantiated in Eqs. (11)–(15) and collectively define the algorithmic novelty of the proposed method.
Algorithmic Comparison Between Standard GP-EHVI and the Proposed Direct GP-EHVI
GP, Gaussian process; GP-EHVI, GP-expected hypervolume improvement; HV, hypervolume.
Figure 3 illustrates the evolution of the Pareto front in the objective space during the MOBO process. The three axes correspond to the three optimization objectives (surface roughness Sa, defect area percentage, and warp RMS). The Pareto fronts generated at each iteration are color-coded with a gradient from light to dark, indicating the progression of optimization. This figure provides an intuitive view of the optimization process: the initial Pareto solutions are widely scattered, and as iterations proceed, the front gradually converges toward the true optimal region, highlights the effectiveness of the multi-objective GP model in guiding the search, improving the distribution of solutions, as well as the trade-offs among different objectives.

Schematic illustration of the Pareto front evolution in the objective space.
Results and Discussion
Comparative algorithm configuration
To ensure a fair evaluation, all three algorithms employ GP regression with the Automatic Relevance Determination squared exponential kernel as the surrogate model, eliminating model capacity as a confounding factor. The GP-based methods (standard GP-EHVI and direct GP-EHVI) operate under an identical budget of 200 evaluations (50 initial Latin Hypercube Sampling samples plus 50 iterations with batch size 3), reflecting the practical constraint of expensive SLM experiments. NSGA-II is allocated 18,000 evaluations to approximate its converged performance under unconstrained conditions, serving as a reference for global search capability rather than a direct competitor under limited budgets.
The critical differences lie in how the surrogate model output is utilized (Table 4). Standard GP-EHVI represents the conventional pipeline in which GP predictions feed passively into the EHVI estimator, with decision-space diversity constraints and a fixed weight (w = 0.5). Direct GP-EHVI establishes three explicit coupling channels: C1 adds an explicit uncertainty bonus (λ = 0.1) to the acquisition function; C2 incorporates GP-predicted objective-space diversity into batch selection; and C3 adaptively regulates the diversity weight according to the HV convergence rate. This comparison isolates the contribution of the proposed coupling architecture from surrogate model quality and evaluation budget effects.
Key Algorithm Parameters for Compared Methods
Comparison of pareto fronts
As illustrated in Figure 4, the Pareto solution distributions obtained by direct GP-EHVI, standard GP-EHVI, and NSGA-II for the SLM process parameter optimization problem are compared from both the perspectives of multi-objective coordination performance and manufacturing-related physical interpretability.

Comparison of Pareto fronts obtained by three optimization methods.
First, from the three-dimensional Pareto front distributions in the objective space (Fig. 4a), it can be observed that the solution set produced by direct GP-EHVI exhibits significantly broader coverage across the three objectives—surface roughness Sa, defect area fraction, and warpage RMS—and approaches the ideal Pareto region more closely along multiple directions. In contrast, the Pareto solutions generated by the standard GP-EHVI method are relatively clustered, mainly concentrated in regions with moderate Sa and moderate warpage RMS, indicating limited exploration capability toward extreme low-defect or low-warpage solutions. The Pareto front obtained by NSGA-II shows both a smaller number of solutions and lower distribution density, with noticeable gaps in regions characterized by high Sa or low defect fraction. These observations suggest that, under a surrogate-model-driven framework, direct GP-EHVI achieves a more effective balance between global exploration and local exploitation, thereby yielding a more representative set of multi-objective trade-off solutions.
Second, the two-dimensional relationship between Sa and defect area fraction (Fig. 4b) reveals that all three methods capture an overall increasing trend of defect fraction with increasing surface roughness. This behavior is consistent with the physical mechanisms of SLM, where higher energy input intensifies melt-pool instability and increases sensitivity to pore formation. However, direct GP-EHVI is capable of identifying solutions with lower defect fractions in the low-Sa region, while also exploring solutions with relatively higher defect fractions at high Sa that may still hold practical relevance within specific process windows. This indicates a higher degree of solution diversity. By comparison, the solutions obtained by standard GP-EHVI are distributed along an approximately single trend line, reflecting a limited ability to characterize nonlinear trade-off relationships between objectives. The NSGA-II solutions are primarily concentrated in intermediate regions, making it difficult to effectively explore potentially optimal solutions associated with extreme process conditions.
Finally, the relationship between Sa and warpage RMS (Fig. 4c) shows a clear decreasing trend of warpage RMS with increasing Sa, suggesting that, within a certain range, higher energy input associated with increased surface roughness contributes to improved melt-pool stability and reduced macroscopic thermal deformation. In this objective plane, direct GP-EHVI not only covers a wider range of Sa values but also maintains favorable control over Sa in the low-warpage region, demonstrating its superiority in addressing SLM quality optimization problems dominated by thermo-mechanical coupling effects. In contrast, both standard GP-EHVI and NSGA-II exhibit insufficient solution continuity and resolution in the low-warpage region.
In summary, direct GP-EHVI demonstrates superior Pareto front approximation accuracy, enhanced solution diversity, and stronger capability in capturing objective trade-offs for multi-objective SLM process parameter optimization. The proposed method more comprehensively reveals the coupled relationships among surface quality, defect formation, and warpage deformation, thereby providing more informative and practically valuable decision support for subsequent process window identification and engineering applications.
Convergence performance analysis
To compare the optimization efficiency and convergence stability of different methods under limited evaluation budgets, the HV indicator was adopted as the primary performance metric. Since the proposed MOBO framework is designed for data-limited optimization problems, rapid improvement of Pareto-front quality with a small number of evaluations was regarded as the main evaluation criterion.
Considering the stochastic nature of BO, both the standard GP-EHVI framework and the proposed improved MOBO method were independently executed using five different random seeds. The mean HV values and corresponding 95% confidence intervals (95% CIs) were calculated to evaluate convergence stability and statistical reliability.
NSGA-II was employed as a representative evolutionary multi-objective optimization baseline, with a population size of 120 and 150 generations, corresponding to ∼18,000 function evaluations. Since evolutionary algorithms generally require significantly larger evaluation budgets than surrogate-assisted BO methods, the HV convergence analysis in Figure 5 mainly focuses on the comparison between the standard GP-EHVI framework and the proposed method, while NSGA-II serves as a reference baseline for global Pareto search capability.

Mean HV convergence curves with 95% CI. CI, confidence interval; HV, hypervolume.
As illustrated in Figure 5, the proposed method consistently achieves higher HV values than the standard GP-EHVI framework throughout the optimization process. In addition, the relatively narrow 95% CIs demonstrate good convergence stability and repeatability across multiple runs. Although the standard GP-EHVI method shows rapid HV improvement during the early stage, its convergence rate gradually decreases afterward. In contrast, the proposed method maintains sustained HV improvement, especially during the mid-to-late optimization stages, indicating a better balance between exploration and exploitation under constrained evaluation budgets.
To further evaluate optimization repeatability, overlaid Pareto fronts from five independent runs are presented in Figure 6 using semi-transparent visualization under identical coordinate axes. The obtained Pareto fronts exhibit relatively concentrated distributions and similar convergence trends, demonstrating good robustness against random initialization effects.

Overlayed Pareto fronts from five runs.
Analysis of solution set distribution characteristics and diversity
In addition to convergence performance, the distribution characteristics of the solution set in the objective space are equally critical to the engineering applicability of multi-objective optimization results. It should be noted that, unlike sequential-sampling-based BO methods, NSGA-II is a population-based evolutionary algorithm in which a large number of individuals are eliminated during generational evolution. As a result, it is difficult to construct a unified solution set that directly corresponds to all evaluated samples, as is typically done in BO. Therefore, for the NSGA-II method, only the final Pareto-optimal solution set is analyzed and compared in this study. Figure 7 compares the distributions of Pareto-optimal solutions obtained by the direct GP-EHVI and the standard GP-EHVI methods in the objective space.

Density distributions of Pareto-optimal solutions obtained by the two methods in the objective space.
Figure 7 presents a comparative visualization of the solution distributions obtained by the standard GP-EHVI method (left column) and the direct GP-EHVI method (right column) for the three-objective optimization problem. The gray scatter points represent all evaluated solutions generated during the optimization process, while the colored scatter points denote the final Pareto-optimal solution set. The three groups of subplots illustrate, respectively, the overall distribution in the three-dimensional objective space and the pairwise objective projections, enabling a comprehensive comparison of the optimization behaviors of the two methods.
First, from the three-dimensional Pareto front distributions in the objective space (Fig. 7a and b), both methods are capable of effectively approximating the multi-objective trade-off surface and forming continuous sets of nondominated solutions. However, compared with the standard GP-EHVI, the Pareto solutions obtained by direct GP-EHVI exhibit a more uniform distribution in the objective space and a smoother, more continuous structure along the front direction. In particular, improved coverage is observed in regions characterized by low defect area percentage and low warp RMS. This indicates that direct GP-EHVI provides greater stability and global consistency in characterizing multi-objective trade-off relationships.
Furthermore, from the two-dimensional projection of Sa versus defect area percentage (Fig. 7c and d), both methods successfully capture the pronounced positive correlation between surface roughness and defect area percentage. In comparison, the Pareto solutions generated by direct GP-EHVI are more compact and exhibit stronger monotonicity in this projection plane, with solution points distributed almost along a smooth curve. This suggests that the method more accurately distinguishes dominated and non-dominated regions during optimization, thereby reducing redundant sampling and local perturbations. In contrast, the solution set obtained by the standard GP-EHVI displays noticeable dispersion in the high-Sa region, indicating that further improvement is possible in exploiting predictive uncertainty under complex nonlinear trade-offs.
Finally, in the Sa–warp RMS projection (Fig. 7e and f), both methods reveal the typical trade-off relationship in which warp RMS decreases as surface roughness increases. Nevertheless, the Pareto solutions obtained by direct GP-EHVI show greater extension in the low-RMS region, continuously exploring improved warp-control solutions. By contrast, the standard GP-EHVI solutions tend to converge prematurely in this region, resulting in a relatively concentrated distribution of solution points. This phenomenon demonstrates that direct GP-EHVI achieves a more effective balance between exploration and exploitation, exhibiting stronger global search capability and effectively avoiding premature aggregation in local optimal regions.
Overall, although the two methods yield broadly consistent Pareto front shapes, direct GP-EHVI demonstrates clear advantages in terms of solution set uniformity, boundary coverage, and stable characterization of multi-objective trade-off relationships. These results indicate that, through the direct coupling and enhancement of the surrogate model and acquisition strategy, direct GP-EHVI can achieve higher-quality Pareto front approximations under a limited evaluation budget, thereby providing more robust decision support for multi-objective optimization of SLM process parameters.
Distribution of process parameters and validation of energy rationality
To further evaluate the differences among the Pareto solutions obtained by different optimization algorithms in terms of physical feasibility and energy input rationality, this section analyzes their characteristics from two perspectives: the distribution in the real process parameter space and the distribution of VED, as shown in Figures 8 and 9.

Comparison of Pareto solution distributions in the real process parameter space.

Boxplot comparison of VED distributions for Pareto solutions. VED, volumetric energy density.
Figure 8 illustrates the distribution of Pareto solutions in the real process parameter space for the three algorithms, together with the two-dimensional LP–SS projection. It can be observed that the Pareto solutions obtained by the standard GP-EHVI model tend to cluster near the boundaries of the design space. This phenomenon mainly arises because, under limited sampling conditions, the surrogate model exhibits relatively high uncertainty at the edges of the design space. Since the EHVI acquisition function accounts for both the predicted mean and variance when calculating the expected improvement, regions with higher uncertainty are preferentially explored, thereby reinforcing the search toward boundary areas. Consequently, this boundary aggregation reflects the sensitivity of the search strategy to uncertain regions rather than the true physical concentration of optimal process parameters.
In contrast, the Pareto solutions obtained by NSGA-II and direct GP-EHVI exhibit more reasonable distributions. Specifically, the NSGA-II solutions are primarily concentrated in regions of relatively low scan speed and low laser power, with a comparatively narrow overall distribution range. This indicates that, under the same evaluation budget, NSGA-II tends to converge to local parameter regions, resulting in limited coverage of the design space. By comparison, the direct GP-EHVI model maintains solutions within a reasonable energy input range while achieving a more balanced and broader parameter distribution, demonstrating enhanced global exploration capability and improved sample utilization efficiency.
Figure 9 further presents the VED distributions of the Pareto solutions obtained by different algorithms. The direct GP-EHVI model exhibits a moderate distribution range, with most solutions concentrated within the engineeringly acceptable energy window. The median remains stable, and no significant extreme values are observed, indicating that the method effectively regulates energy input while achieving performance optimization, thereby balancing densification and the risk of overheating.
In contrast, the VED distribution of NSGA-II is relatively concentrated but narrower, suggesting that its search results are largely confined to empirically stable regions. Although this conservative tendency contributes to process stability, it may restrict the potential for further performance improvement.
Due to the boundary aggregation observed in the parameter space, the standard GP-EHVI model shows a wider VED spread, with some solutions exceeding the reasonable energy window. This reflects an over-exploration of high-uncertainty regions, leading to deviations from practically stable processing conditions.
Overall, the direct GP-EHVI model achieves a better balance between parameter distribution and energy rationality, demonstrating superior engineering applicability.
To elucidate the interaction effects among LP, SS, and HS, Figure 10 presents the response surface mapping between process parameters and surface quality metrics. Surface roughness increases with decreasing LP and increasing SS, as reduced energy input lowers melt pool temperature and prevents complete powder melting. Roughness is particularly sensitive to HS variation, where large spacing reduces track overlap and produces stripe patterns. The highest roughness occurs when low LP coincides with large HS, where insufficient energy and poor overlap jointly degrade melt pool stability. Defect percentage follows similar trends, with inadequate melting elevating defects under low energy input. Warp deformation is maximized by high LP with low SS through intensified thermal gradients, yet becomes scan-speed dominated at higher speeds as reduced per-length heat input mitigates thermal stress. Notably, warp distributes bimodally, with moderate HS and low SS yielding high warp while parameter extremes produce lower values, indicating coupled energy-overlap effects.

Multi-objective response surface mapping between SLM process parameters and surface quality indicators. SLM, selective laser melting.
A critical finding is that similar VED values can arise from markedly different parameter combinations with divergent quality outcomes. For example, two solutions at approximately 50 J/mm³—(190 W, 1000 mm/s, 95 μm) versus (150 W, 650 mm/s, 110 μm)—exhibit opposite trade-offs: the former achieves lower Sa and defects but higher warp through a larger, stable melt pool with intensified thermal gradients, while the latter shows higher Sa and defects but lower warp through a smaller, less stable melt pool with reduced thermal stress. This demonstrates that VED alone is insufficient to predict surface quality; the specific parameter combination determines spatiotemporal energy distribution and consequent melt pool dynamics.
Figure 11 further confirms this relationship through the VED scatter distributions. Although Sa and defect percentage generally decrease while warp increases with increasing VED, substantial scatter remains at equivalent VED levels, indicating that surface quality is not determined solely by energy density. Instead, different parameter combinations can produce significantly different melting behavior, melt pool stability, and thermal accumulation characteristics even under similar VED conditions, ultimately leading to variations in surface morphology and overall surface quality.

Distribution of surface quality indicators as a function of VED.
Experimental validation
To evaluate the feasibility and effectiveness of the multi-objective optimization results in practical manufacturing, experimental validation was conducted based on 10 representative optimal process parameter sets selected from the Pareto front. As shown in Table 5, all selected Pareto-optimal solutions achieved stable fabrication without evident process failure or extreme anomalies, indicating that the optimized parameter combinations fall within a manufacturable processing window.
Pareto-Optimal Process Parameter Combinations and Their Corresponding Surface Quality Metrics
From the perspective of performance metrics, the experimental results exhibit consistent variation trends among surface roughness, defect ratio, and warpage, which agree well with the optimization outcomes. Specifically, in the low-roughness region (Samples 1–4), the measured Sa values are maintained within ∼6–7 μm, while the defect ratio and warpage remain at moderate levels. In contrast, as the process parameters shift toward higher scanning speeds or lower energy input conditions (Samples 8–10), both surface roughness and defect ratio increase significantly, whereas warpage tends to decrease. This behavior is consistent with the trade-off relationships reflected by the Pareto front, demonstrating that the optimization results can effectively capture the intrinsic coupling mechanisms among different performance objectives.
The quantitative comparison in Table 6 reveals that the prediction accuracy varies significantly across objectives and samples, and these variations must be interpreted in the context of both the physical measurement process and the surrogate model characteristics.
Comparison Between Predicted and Experimental Values of Optimal Process Parameter Combinations
APE, absolute percentage error; |Pred. − Exp.|/Exp. × 100%.
For Sa, the absolute percentage error (APE) values are generally low (<10%) for most samples, with the exception of Sample 6 (APE = 13.2%). The small deviations observed for Samples 1–5 (APE < 8%) indicate that the GP surrogate effectively captures the surface morphology response within the interpolation-dominant region of the parameter space. The elevated error for Sample 6 (predicted 7.75 μm vs. experimental 8.77 μm) corresponds to the highest laser power (200 W) and scan speed (1197 mm/s) in the validation set, representing an extreme operating condition where melt-pool dynamics deviate from the training data distribution.
For defect percentage, the predictions exhibit substantial relative errors in the low-defect regime (Samples 1–4, APE = 37.6–51.7%). This phenomenon is physically attributable to the inherent difficulty in predicting rare defects: When the true defect percentage falls below 2%, even small absolute deviations (e.g., 1.27 percentage points for Sample 1) translate into large percentage errors. Importantly, the predicted rankings remain correct—Samples 1–4 are consistently predicted and experimentally confirmed as the low-defect group—indicating that the model captures the qualitative discrimination capability despite quantitative imprecision in the extreme low-defect regime. For Samples 6–10 where defect percentages exceed 4%, the APE drops to below 15%, confirming improved prediction reliability in the moderate-to-high defect range.
For warp RMS, the predictions are generally accurate (APE < 10%) for 8 of the 10 samples, with two notable exceptions: Sample 6 (APE = 34.3%, predicted 9.66 μm vs. experimental 12.97 μm) and Sample 8 (APE = 10.9%, predicted 11.21 μm vs. experimental 10.11 μm). The substantial underestimation for Sample 6 warrants specific discussion. This sample corresponds to the highest energy density (VED = 62.1 J/mm³) and the most extreme thermal conditions in the validation set. Under such high-power, high-speed conditions, the thermal gradient and residual stress distribution deviate significantly from the interpolation-dominant region of the training data, leading to extrapolation uncertainty in the GP surrogate. The physical mechanism involves accelerated cooling rates and asymmetrical thermal contraction that amplify warpage beyond the surrogate’s predictive capacity. This observation highlights a fundamental limitation of surrogate-based optimization: Predictions at parameter space boundaries carry higher uncertainty, and experimental validation is essential to anchor these extreme predictions.
Overall, the experimental validation confirms three critical aspects of the optimization framework. First, despite local quantitative deviations—particularly for boundary samples (Sample 6) and low-defect predictions (Samples 1–4)—the predicted Pareto-optimal solutions consistently fall within the manufacturable process window and exhibit the expected trade-off behavior. Second, no “predicted optimal but experimentally infeasible” solutions were encountered, confirming the physical plausibility of the Pareto front. Third, the model successfully identifies the correct qualitative rankings across all three objectives, which is the primary requirement for multi-objective optimization guidance. The identified large errors are not dismissed as minor deviations but are explicitly attributed to identifiable physical and modeling factors: boundary extrapolation effects for extreme parameters and signal-to-noise limitations for rare-defect prediction. These insights inform future work on adaptive sampling strategies targeted at high-uncertainty regions, such as incorporating expected improvement of prediction variance criteria to selectively sample boundary regions during optimization.
It is important to interpret the observed deviations in the context of the inherent variability of the SLM process. Previous studies have shown that even under identical processing conditions, measurable variations in surface morphology and defect characteristics may still arise due to the stochastic behavior of melt-pool evolution, powder spreading, and thermal fluctuations.36,37 Consequently, a certain degree of deviation between predicted and experimentally measured values is unavoidable in practical LPBF/SLM systems. In this study, the validation results are therefore interpreted primarily from the perspective of qualitative trend consistency and engineering feasibility rather than strict point-wise statistical accuracy. Overall, the observed prediction deviations remain broadly consistent with the variability characteristics commonly reported in the literature, indicating that the proposed optimization framework still possesses meaningful engineering reference value despite unavoidable process uncertainty.
Conclusion
This study proposed a multi-objective BO approach based on direct GP-EHVI for the collaborative optimization of SLM process parameters. By introducing a direct coupling improvement between the surrogate model and the acquisition function, the proposed method enables a refined characterization of the coupled relationships among key quality indicators, including surface roughness, defect formation, and warpage deformation. Under a limited function evaluation budget, the method demonstrates superior Pareto front approximation capability, higher solution diversity, and more stable HV improvement efficiency compared with conventional GP-EHVI and NSGA-II, particularly maintaining continuous advancement of the Pareto front during the middle and later optimization stages.
The results indicate that the proposed method exhibits clear advantages in solution uniformity, boundary coverage, and stable representation of multi-objective trade-offs. It achieves a favorable balance between rational parameter space distribution and controllable energy input. Moreover, the framework significantly reduces reliance on large-scale experimental campaigns and computationally expensive numerical simulations, enabling the identification of practically relevant process parameter combinations within a limited number of trials.
Overall, the proposed approach outperforms the benchmark methods in both optimization performance and engineering applicability, providing a robust and scalable technical pathway for intelligent multi-objective optimization of SLM process parameters. It offers valuable support for efficient decision-making and quality control in complex additive manufacturing processes.
Future research will extend the proposed method to higher-dimensional parameter spaces and more complex processing conditions. Physics-informed modeling strategies will be incorporated to enhance robustness and generalization under small-sample scenarios. In addition, the integration of in-situ monitoring data and feedback control mechanisms will be explored to promote the development of adaptive and closed-loop intelligent optimization for SLM processes.
Authors’ Contributions
All authors contributed to the conception and design of the study. Material preparation, data collection, and analysis were performed by S.C. and H.C. The experiment was carried out by S.C. and Y.Z. The original article was organized and written by H.C., S.C., and Y.Z. S.D. provided methodology and supervised the completion of the article. All authors read and approved the final article.
Footnotes
Author Disclosure Statement
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this article.
Funding Information
This work was financially supported by the grants from the
