Integrated Taguchi–gray relational analysis for multi-objective optimization of fused deposition modeling process parameters to balance carbon emissions and ultimate tensile strength

Abstract

This study optimizes fused deposition modeling (FDM) parameters to resolve the trade-off between carbon emission reduction and ultimate tensile strength (UTS). Using polylactic acid (PLA), a Taguchi $L_{16}$ orthogonal array was employed to evaluate five factors: layer height, nozzle temperature, bed temperature, printing speed, and infill rate. Analysis of variance (ANOVA) and response surface methodology (RSM) identified infill rate and layer height as the most significant factors. To address the inherent conflict between environmental impact and mechanical performance, gray relational analysis (GRA) was implemented for multi-objective optimization. The optimal configuration (0.15 mm layer height, 200 °C nozzle temperature, 70 °C bed temperature, 50 mm/s speed, and 30% infill rate) increased UTS by 18.69% relative to the minimum-emission setting, and reduced carbon emissions by 15.88% compared to the maximum-strength setting. This framework provides a critical reference for achieving low-carbon, high-efficiency additive manufacturing in industrial applications.

Introduction

With the global emphasis on sustainable manufacturing and carbon neutrality goals, reducing energy consumption and carbon emissions during manufacturing has become a core issue in both industrial and academic research. ¹ In practical applications, identifying high-emission stages through scientific carbon footprint auditing, coupled with real-time energy efficiency monitoring via smart meters, can assist enterprises in optimizing decarbonization strategies and strengthening their competitive position. ²

Additive manufacturing (AM) is widely regarded as a technology with environmental and efficiency advantages over certain conventional processes.^3,4 However, its actual environmental performance depends heavily on specific process modalities and parameter configurations; without proper control, AM may not perform as anticipated in terms of energy consumption or carbon footprint. ⁵

Building upon the aforementioned literature, this study utilizes a Taguchi $L_{16}$ orthogonal array in conjunction with real-time machine energy monitoring and carbon emission conversion. The goal is to systematically examine the effects of four primary process parameters—printing speed, layer thickness, infill rate, and printing temperature—on the carbon footprint and tensile strength of individual components. This research extends the application of the Taguchi method to contexts where environmental performance and material properties are equally prioritized.

Theoretical background

The quantification of energy expenditure and the resulting carbon footprint in additive manufacturing (AM) has evolved from static G-code-based prediction models to dynamic, multi-dimensional life cycle assessment (LCA) frameworks.^6,7 Recent literature emphasizes that the environmental performance of fused deposition modeling (FDM) is highly sensitive to real-time power fluctuations during the thermal extrusion phase. ⁸ Emerging studies from 2025 underscore that integrating high-precision real-time monitoring with LCA metrics provides a more granular assessment of AM sustainability, as energy consumption during the polymer melting stage can account for the majority of a part’s total carbon emissions.^9–11 A comprehensive overview of these processes indicates that while AM offers transformative potential, the industry faces persistent challenges in stabilizing material-dependent energy footprints. ¹²

Crucially, the energy demand during this thermal phase is inextricably linked to the polymer’s rheological behavior and the material’s thermal processing window.^13,14 Factors such as nozzle temperature and layer height directly govern the melt viscosity and the cumulative thermal history of the specimen, which dictate the degree of molecular diffusion and interlayer fusion across deposition interfaces. This relationship becomes even more complex with the development of novel PLA-based biocomposites, such as those reinforced with Henna ¹⁵ or Wood,^16,17 where the introduction of organic fillers significantly alters the thermal processing requirements and the evolution of mechanical properties like flexural strength. ¹⁸ Ensuring sufficient polymer chain mobility without inducing thermal degradation is essential for maximizing ultimate tensile strength (UTS), while simultaneously managing the heating duty cycle to minimize the cumulative environmental impact. ¹⁹

The Taguchi method is recognized for its proficiency in analyzing multi-factor interactions and enhancing process robustness through signal-to-noise (S/N) ratio analysis.^20,21 While its efficacy has been validated in diverse manufacturing domains—such as automotive forging design ²² and the reduction of defects in family mold production ²³ —traditional applications have primarily focused on mechanical or dimensional quality characteristics. Extending this methodology to prioritize environmental sustainability metrics, specifically carbon emissions, as a response variable remains an active research frontier that requires further investigation into its multi-criteria viability.^24–26 Consequently, the integration of computational intelligence—such as ANOVA combined with backpropagation neural networks ²⁷ or supervised learning ¹⁶ —is increasingly required to model the complex mechanical behavior and environmental trade-offs of advanced composite materials. This hybrid approach allows for a more granular understanding of how varied process configurations influence both structural integrity and the instantaneous energy shifts within the manufacturing cycle.

A fundamental challenge in FDM optimization lies in the inherent conflict between mechanical performance and environmental objectives. Recent research by Salehi and Movahhedy ²⁸ highlighted these antagonistic characteristics, demonstrating that parameters optimized for structural integrity, such as elastic modulus, often compromise other quality indicators like surface topography. To resolve these Pareto-optimal trade-offs, the hybridization of response surface methodology (RSM) and gray relational analysis (GRA) has emerged as a high-fidelity solution for reconciling quality and efficiency in advanced manufacturing.^29,30 This multi-criteria approach is particularly instrumental when balancing complex variables, such as reinforcing phases in composite materials with surface integrity.^31,32

Research method

This study integrates the Taguchi method, gray relational analysis (GRA), and response surface methodology (RSM) to overcome the inherent limitations of using these techniques independently. While the Taguchi method efficiently minimizes experimental runs, it is confined to discrete, single-objective evaluations. GRA successfully resolves multi-objective conflicts between ultimate tensile strength and carbon emissions, but it cannot predict outcomes outside the predefined experimental grid. By applying RSM to the Taguchi-GRA derived gray relational grade (GRG), this tripartite framework transcends discrete boundaries to generate a continuous polynomial prediction model. This synergistic approach significantly enhances optimization reliability, enabling the precise identification of the global Pareto-optimal parameter configuration while avoiding the experimentally prohibitive costs of a full-factorial RSM design.

This study is organized into four chapters, and the corresponding research framework is illustrated in Figure 1.

Figure 1.

Research flowchart.

Research model

The ASTM D638 standard is widely used to evaluate the tensile properties of various materials. Accordingly, tensile testing in this study was conducted in accordance with the ASTM D638-10 Type I specification. ³³ The detailed dimensions and specifications of the specimens are presented in Figure 2.

Figure 2.

Research model. ¹⁹

Experimental equipment

Specimens were fabricated using a DUAL300 FDM system (Linchen Industrial Co., Ltd., Taiwan), as shown in Figure 3. Polylactic acid (PLA) was selected as the printing material. The system integrates a precision extrusion mechanism, a dual-zone thermal control unit (nozzle and build plate), and an active cooling system to ensure dimensional stability during rapid solidification. The fabrication process followed G-code trajectories that defined the layer-by-layer deposition paths.

Figure 3.

3D printer. ¹⁹

Polylactic acid (PLA) was employed as the primary printing material to fabricate the experimental models and specimens. As a bio-based polymer, PLA’s intrinsic properties—including hygroscopicity, thermal stability, and morphological integrity—are highly sensitive to environmental fluctuations. These sensitivities may introduce inconsistencies during extrusion, thereby affecting the dimensional accuracy and impact toughness of the fabricated components. Existing research on PLA indicates that a finer layer height enhances tensile strength and that internal infill patterns provide superior structural support. ³⁴ Consequently, the PLA filament used in this research, as shown in Figure 4, was conditioned and processed under controlled environmental parameters to ensure data reliability and structural consistency.

Figure 4.

3D printing filament. ¹⁹

The slicing and process preparation were performed using PING Slicer, a proprietary software suite developed by the printer manufacturer. This platform facilitated the configuration of process parameters and the translation of 3D CAD models into G-code instructions, as illustrated by the generated toolpaths in Figures 5 and 6. Given that process variables directly influence the dimensional fidelity, surface topography, and mechanical performance of the printed specimens, precise control within the software environment is essential to ensure experimental reproducibility.

Figure 5.

GUI for process parameter configuration.

Figure 6.

Model diagram in the software.

The fabricated standard specimens were subjected to uniaxial tensile testing using a universal testing machine (UTM), as illustrated in Figure 7. Upon securing the specimens within the grips, a tensile load was applied at a constant crosshead speed of 5 mm/min until fracture occurred. This testing protocol enabled the characterization of key mechanical properties, including the yield strength (σ_y), ultimate tensile strength (UTS), and ductility (percentage elongation).

Figure 7.

Universal testing machine. ¹⁹

Energy monitoring system

To quantify the environmental footprint, a real-time energy monitoring system was implemented, consisting of three primary modules:

• 1. Power sensing module:

  A high-precision single-phase power analyzer (accuracy of ±0.1) was integrated into the equipment’s power input. This module facilitates synchronous measurement of real-time voltage, current, and active power, enabling precise recording of instantaneous power fluctuations and cumulative energy consumption.

• 2. Data acquisition and logging system:

  The data acquisition (DAQ) platform, comprising integrated hardware-software interfaces, was designed to monitor and record total energy consumption (kWh) throughout the printing cycle, as illustrated in Figure 8. This configuration ensured high-fidelity data collection for subsequent energy efficiency analysis.

• 3. Carbon emission calculation module:

  Following the principles of life cycle assessment (LCA), this study established a quantifiable carbon emission model ⁵ by integrating material throughput, equipment energy usage, and the regional electricity carbon emission factor for Taiwan. ³⁵ To ensure a balance between analytical rigor and practical operability, the system boundary focused on two primary emission sources: the material phase and the energy phase.

Figure 8.

Energy consumption system. ³⁶ (The red square box highlights the real-time active power measurement value of 160.0 W during the printing phase).

To ensure the fidelity and validity of this LCA framework, the carbon emission model must be grounded in precise, boundary-specific emission factors. For the material phase, carbon emissions are determined by the product of total material consumption and its unit emission factor. Specifically, the carbon emission factor for polylactic acid (PLA) was adopted from the comprehensive LCA inventory established by Vink and Davies. ¹³ This factor accounts for cradle-to-gate environmental impacts, encompassing agricultural feedstock production, fermentation, and polymerization processes, thereby ensuring alignment with international biopolymer assessment standards.

In the energy phase, emissions were quantified by monitoring real-time power demand during the FDM process. Recognizing that the electricity emission factor is highly sensitive to the regional power generation mix, this study avoided using generic global averages. Instead, cumulative energy consumption was converted using the official regional electricity carbon emission factor published by the Bureau of Energy, Ministry of Economic Affairs (Taiwan Power Company). ³⁵ By explicitly integrating peer-reviewed material parameters with specific regional grid intensities, the accuracy and spatial-temporal relevance of the proposed environmental evaluation model are rigorously validated. The total carbon emissions ( $C_{\mathrm{total}}$ ) were calculated using equations (1) through (3):

$$C_{\mathrm{total}}=C_{\mathrm{materia}}+C_{\mathrm{energy}}$$

(1)

$$C_{\mathrm{material}}=m\times {\mathrm{EF}}_{\mathrm{material}}$$

(2)

$$C_{\mathrm{energy}}=P\times t\times {\mathrm{EF}}_{\mathrm{electricity}}$$

(3)

Where:

 m is the actual mass of the material consumed (kg);

  ${\mathrm{EF}}_{\mathrm{material}}$ is the cradle-to-gate carbon emission factor per unit of material (kg  ${\mathrm{CO}}_2$ e/kWh);

 P represents the average power demand of the 3D printer during the printing process (kW);

 t is the total printing duration (h);

  ${\mathrm{EF}}_{\mathrm{electricity}}$ denotes the electricity carbon emission factor for Taiwan (kg  ${\mathrm{CO}}_2$ e/kWh).

According to the Administration of Energy, Ministry of Economic Affairs (2025), the electricity carbon emission factor for 2024 is 0.474 kg  ${\mathrm{CO}}_2$ e/kWh. ³⁵ For the material phase, the emission factor for PLA was set at 1.80 kg  ${\mathrm{CO}}_2$ e/kWh, derived from established life cycle assessment (LCA) databases for bio-based polymers. ³² These coefficients were integrated into a quantitative model to evaluate the impact of process configurations, enabling the identification of parameters that reconcile low-carbon production with enhanced mechanical performance.

Results and discussion

This section presents the evaluation of carbon emissions during the fabrication process and the tensile strength of the 3D-printed specimens. To ensure statistical significance and minimize experimental error, each experimental group was subjected to four independent replicates. The mean value ( $y_i$ ) derived from these four measurements was calculated and used as the representative data point for each configuration. Subsequently, these averaged results were processed and analyzed using the Taguchi method to quantify the influence of various processing parameters on the defined quality characteristics.

Taguchi method

The Taguchi method provides a robust experimental framework by utilizing orthogonal arrays (OAs) to systematically organize control factors and their respective levels. This methodology integrates S/N ratio analysis and ANOVA to evaluate the sensitivity of processing parameters to the defined quality characteristics. By effectively minimizing the influence of noise factors and experimental fluctuations—which were reduced in this study by calculating mean values from four independent replicates—this framework facilitates the identification of an optimal parameter configuration with a significantly reduced number of experimental trials.

Taguchi method of ultimate carbon emissions

In this study, five critical control factors governing the quality characteristics were identified, each encompassing four distinct levels. The selection of these independent variables for the Taguchi array was governed by a strict multi-objective criterion, grounded in a comprehensive review of established literature^37–40: the parameters must significantly influence both structural integrity and environmental impact. Consequently, variables like raster angle were excluded; although it affects directional strength, it demonstrates negligible impact on material consumption and process duration, contributing no meaningful variance to carbon emissions. Similarly, cooling rate was designated as a constant control, as maximal cooling is a thermodynamic prerequisite for PLA to prevent thermal warping that would invalidate mechanical tests. Therefore, layer height, nozzle temperature, bed temperature, printing speed, and infill rate were isolated as the core thermal and kinetic drivers.

The specific levels for each parameter were defined based on machine capabilities and material properties, as detailed in Table 1. These discrete levels were strategically established to encompass the viable processing window of polylactic acid (PLA). The thermal boundaries—nozzle temperature (190 °C–220 °C) and bed temperature (50 °C–80 °C)—were calibrated to optimize melt viscosity and polymer chain mobility. This range ensures sufficient molecular diffusion for interlayer fusion while strictly avoiding thermal degradation and excessive energy expenditure. Furthermore, geometric and kinetic factors such as layer height (0.10–0.25 mm) and print speed (40–70 mm/s) were selected to modulate the cumulative thermal history and cooling gradients of the specimen. Finer layer heights and lower speeds enhance interfacial bonding (benefiting UTS) but significantly prolong the duty cycle of the machine’s heating elements, thereby inflating the energy-phase carbon footprint. Finally, infill density (30%–90%) was defined to capture the linear scalability of structural robustness against material volume.

Table 1.

Control factor level.

Control factor	Level 1	Level 2	Level 3	Level 4
A. Layer height (mm)	0.10	0.15	0.20	0.25
B. Nozzle temperature (°C)	190	200	210	220
C. Bed temperature (°C)	50	60	70	80
D. Print speed (mm/s)	40	50	60	70
E. Infill density (%)	30	50	70	90

Together, these parameters and their defined levels construct a robust design space. They directly dictate material throughput and heating duty cycles, providing the essential variance to effectively model the Pareto-optimal trade-offs between mechanical integrity and environmental sustainability without inducing macroscopic printing failures.

Post-fabrication, carbon emissions were quantified based on cumulative power consumption, while the specimens were subjected to tensile testing to assess mechanical integrity. For the Taguchi analysis, carbon emissions were evaluated using the “Smaller-the-Better” (STB) criterion to minimize environmental impact, as expressed in equation (4). Conversely, tensile strength was analyzed using the “Larger-the-Better” (LTB) criterion to maximize structural performance, as defined in equation (5). The consolidated experimental results are presented in Table 2.

$$S/N_{\mathrm{SB}}=-10\log \left[\frac{1}{n}\sum _{i=1}^n{y}_i^2\right]$$

(4)

$$S/N_{\mathrm{lB}}=-10\log \left[\frac{1}{n}\sum _{i=1}^n\frac{1}{y_i^2}\right]$$

(5)

where n denotes the total number of experimental trials, and $y_i$ represents the average data for each experimental group.

Table 2.

S/N ratios of carbon emissions for each group.

Exp.	A	B	C	D	E	Carbon_y1	Carbon_y2	Carbon_y3	Carbon_y4	Carbon emissions	$S/{N\prime }_{\mathrm{ce}}$
1	1	1	1	1	1	24.9	24.8	25.9	23.2	24.7	−27.9
2	1	2	2	2	2	35.6	32.4	34.9	34.7	34.4	−30.7
3	1	3	3	3	3	43.6	43.7	42.6	45.3	43.8	−32.8
4	1	4	4	4	4	52.2	51.0	51.3	49.5	51.0	−34.2
5	2	1	2	3	4	39.0	39.4	38.5	37.9	38.7	−31.8
6	2	2	1	4	3	35.8	35.9	36.3	35.2	35.8	−31.1
7	2	3	4	1	2	34.6	34.7	34.1	35.4	34.7	−30.8
8	2	4	3	2	1	32.2	33.6	34.1	31.3	32.8	−30.3
9	3	1	3	4	2	27.7	28.3	27.0	29.0	28.0	−28.9
10	3	2	4	3	1	25.1	25.0	26.6	24.9	25.4	−28.1
11	3	3	1	2	4	38.1	36.7	36.1	40.3	37.8	−31.5
12	3	4	2	1	3	38.0	38.9	39.1	39.2	38.8	−31.8
13	4	1	4	2	3	27.9	25.8	27.8	27.7	27.3	−28.7
14	4	2	3	1	4	34.6	35.0	35.1	33.7	34.6	−30.8
15	4	3	2	4	1	24.9	23.9	23.1	20.1	23.0	−27.2
16	4	4	1	3	2	24.5	23.4	23.4	25.5	24.2	−27.7
Average										33.4	−30.3

As evidenced by the response table (Table 3) and the main effects plot for S/N ratios (Figure 9), a higher S/N ratio indicates reduced variability relative to the target value, effectively correlating to a reduction in carbon emissions. The optimal process configuration for minimizing the carbon footprint was identified as follows: layer height of 0.25 mm, nozzle temperature of 190 °C, bed temperature of 50 °C, print speed of 60 mm/s, and infill density of 30%.

Table 3.

Response table of factors to S/N ratio (carbon emissions).

Control factor	A	B	C	D	E
Level 1	−31.39	−29.32	−29.54	−30.3	−28.38
Level 2	−30.99	−30.17	−30.37	−30.33	−29.54
Level 3	−30.09	−30.61	−30.72	−30.09	−31.1
Level 4	−28.6	−30.98	−30.44	−30.35	−32.06
Effect	2.79	1.66	1.18	0.26	3.68
Rank	2	3	4	5	1
Optimal combination	A4	B1	C1	D3	E1

Figure 9.

Response diagram of factors to the S/N ratio (carbon emissions).

Based on the Delta statistics derived from the factor effects, infill density emerged as the most dominant determinant of the environmental response, whereas print speed demonstrated a comparatively limited influence within the investigated parameter range.

Taguchi method of ultimate tensile strength

The experimental design incorporated five control factors, each with four levels, and detailed parameter configurations are listed in Table 1. This subsection analyzes the tensile strength results obtained from the UTM. Given that maximizing structural integrity was the primary objective, the “Larger-the-Better” (LTB) criterion of the Taguchi method—as defined in equation (5)—was applied for the S/N ratio optimization. The corresponding experimental data and S/N ratios are summarized in Table 4.

Table 4.

S/N ratios of UTS for each group.

Exp.	A	B	C	D	E	UTS_y1	UTS_y2	UTS_y3	UTS_y4	Ultimate tensile strength (MPa)	$S/{N\prime }_{\mathrm{UTS}}$
1	1	1	1	1	1	32.9	32.8	33.9	31.2	32.7	30.3
2	1	2	2	2	2	39.6	37.4	38.9	38.7	38.7	31.7
3	1	3	3	3	3	42.6	42.7	41.6	43.3	42.6	32.6
4	1	4	4	4	4	44.7	43.5	43.8	42	43.5	32.8
5	2	1	2	3	4	40	40.4	39.5	38.9	39.7	32.0
6	2	2	1	4	3	38.3	38.4	38.8	37.7	38.3	31.7
7	2	3	4	1	2	36.6	36.7	36.1	37.4	36.7	31.3
8	2	4	3	2	1	41.2	42.6	43.1	40.3	41.8	32.4
9	3	1	3	4	2	33.2	33.8	32.5	34.5	33.5	30.5
10	3	2	4	3	1	30.1	30	31.6	29.9	30.4	29.7
11	3	3	1	2	4	42.1	40.7	40.1	43.3	41.6	32.4
12	3	4	2	1	3	39	39.9	40.1	40.2	39.8	32.0
13	4	1	4	2	3	35.9	33.8	35.8	35.7	35.3	31.0
14	4	2	3	1	4	38.6	39	39.1	37.7	38.6	31.7
15	4	3	2	4	1	33.4	32.4	31.6	28.6	31.5	30.0
16	4	4	1	3	2	33.5	32.4	32.4	34.5	33.2	30.4
Average										37.4	31.4

According to the analysis of Table 5 and Figure 10, a higher S/N ratio in the “Larger-the-Better” (LTB) characteristic indicates reduced variability and enhanced ultimate tensile strength. The optimal parameter combination was identified as follows: layer height of 0.10 mm, nozzle temperature of 220 °C, bed temperature of 70 °C, print speed of 50 mm/s, and infill density of 90%. Influence analysis indicated that infill density was the most influential factor, whereas layer height and nozzle temperature demonstrated comparatively lower effects. Bed temperature was identified as the least significant factor within the evaluated range.

Table 5.

Response table of factors to S/N ratio (UTS).

Control factor	A	B	C	D	E
Level 1	31.84	30.93	31.18	31.33	30.58
Level 2	31.84	31.2	31.42	31.87	30.99
Level 3	31.13	31.55	31.81	31.16	31.8
Level 4	30.77	31.9	31.17	31.23	32.21
Effect	1.07	0.97	0.64	0.71	1.63
Rank	2	3	5	4	1
Optimal combination	A1	B4	C3	D2	E4

Figure 10.

Response diagram of factors to S/N ratio (UTS).

ANOVA

ANOVA was used to evaluate the sources of error in the experimental dataset and to determine the statistical significance of each control factor using the F-test. By assessing the contribution of individual parameters, the primary determinants of both carbon emissions and tensile strength can be systematically identified.

ANOVA of carbon emissions

To assess the statistical significance of the investigated parameters, ANOVA was performed on the carbon emission data from Section 3.1.1. Utilizing Minitab, the analysis combined the experimental levels from Table 1 with the corresponding performance metrics in Table 2. The resulting F-values and p values characterizing the influence of each control factor on carbon emissions are presented in Tables 6 and 7.

Table 6.

The first S/N ratio variation analysis (carbon emissions).

Factor	SS	DOF	Variance
A. Layer height (mm)	18.3	3	6.1
B. Nozzle temperature (°C)	6.1	3	2.0
C. Bed temperature (°C)	3.1	3	1.0
D. Print speed (mm/s)	0.2	3	0.1
E. Infill density (%)	32.1	3	10.7
Other	*	0	*
Total	59.8	15

*

Pooled into error.

Table 7.

The second S/N ratio variation analysis (carbon emissions).

Second S/N ratio variation analysis
Factor	SS	DOF	Variance	F Statistic (F)	ρ (%)	Confidence	Significant
A	18.3	3	6.1	103.6	30.60%	99.8%	Yes
B	6.1	3	2.0	34.6	10.20%	99.2%	Yes
C	3.1	3	1.0	17.5	5.18%	97.9%	Yes
D	Pooled
E	32.1	3	10.7	181.3	53.68%	99.9%	Yes
Error	0.2	3	0.1	Note: At least 95% confidence
Total	59.8	15

Based on the initial ANOVA results in Table 6 print speed (D) exhibited the least variance and was subsequently pooled into the experimental error to improve the statistical power of the second ANOVA, as presented in Table 7.

Mechanistically, the predominant influence of infill density (53.68%) and layer height (30.60%) on the environmental response can be attributed to their direct governance over the material and energy vectors of the FDM process. Carbon emissions in this framework are bipartite, originating from both material consumption and electricity usage.

Infill density acts as the primary driver for material-phase emissions; higher density configurations exponentially increase the required volume of PLA. Concurrently, a denser internal geometry dictates a more extensive and complex toolpath, thereby prolonging the active extrusion time per layer.

Conversely, layer height serves as a critical geometric constraint that inversely dictates the total process duration. A reduction in layer thickness proportionally increases the total number of discrete layers required to fabricate the component. Because the continuous thermal maintenance of the extrusion nozzle and the build bed dominates the power profile of the FDM system, printing time effectively acts as a multiplier for energy consumption. Consequently, finer layer heights drastically extend the thermal duty cycle of the equipment, cementing layer height and infill density as the most highly sensitive parameters in the life cycle carbon footprint of the printed parts.

The refined analysis indicated that infill density (E) was the most dominant determinant of carbon emissions, accounting for 53.68% of the total percentage contribution. This was followed by layer height (A) at 30.60%, while bed temperature was identified as the least influential factor, contributing only 5.18%. Notably, all investigated factors achieved statistical significance at the 95% confidence level (p < 0.05), confirming that each parameter exerted a statistically significant impact on the environmental footprint of the FDM process.

ANOVA of UTS

This subsection provides a detailed ANOVA of the UTS data derived from the Taguchi experimental framework discussed in Section 3.1.2. The statistical analysis was performed using Minitab by correlating the control factor levels defined in Table 1 with the mean UTS values and their corresponding S/N ratios presented in Table 4.

The primary objective of this analysis was to decompose the total experimental variance into specific sources, thereby ascertaining the statistical significance and percentage contribution of each printing parameter to the mechanical integrity of the specimens. The resulting ANOVA statistics, including F-values and p values, are summarized in Table 8 (initial analysis) and Table 9 (refined analysis).

Table 8.

The first S/N ratio variation analysis (UTS).

Factor	SS	DOF	Variance
A. Layer height (mm)	3.431	3	1.1437
B. Nozzle temperature (°C)	2.147	3	0.7155
C. Bed temperature (°C)	1.059	3	0.353
D. Print speed (mm/s)	1.256	3	0.4188
E. Infill density (%)	6.591	3	2.1969
Other	*	0	*
Total	14.484	15

*

Pooled into error.

Table 9.

The second S/N ratio variation analysis (UTS).

Second S/N ratio variation analysis
Factor	SS	DOF	Variance	F Statistic (F)	ρ (%)	Confidence	Significant
A	3.431	3	1.1437	3.24	23.69%	82.00%	Yes
B	2.147	3	0.7155	2.03	14.82%	71.20%	No
C	Pooled
D	1.256	3	0.4188	1.19	8.67%	55.40%	No
E	6.591	3	2.1969	6.22	45.51%	91.60%	Yes
Error	1.059	3	0.353	Note: At least 80% confidence
Total	14.484	15

Based on the initial ANOVA results presented in Table 8, bed temperature (C) exhibited the smallest variance among the factors. To improve the statistical power for the dominant parameters, bed temperature was pooled into the experimental error term for a subsequent, refined ANOVA, as detailed in Table 7.

As indicated by the refined analysis (Table 9), infill density (E) emerged as the most dominant determinant of mechanical performance, accounting for 45.51% of the total percentage contribution. This was followed by layer height (A) at 23.69% and nozzle temperature (B) at 14.82%, while print speed (D) was found to be the least influential factor, contributing only 8.67%. Notably, except for nozzle temperature and print speed, all other factors achieved a confidence level exceeding 80%, demonstrating a statistically significant impact on the UTS of the 3D-printed specimens within the investigated parameter range.

RSM

RSM is a statistical and mathematical framework that integrates designs of experiments (DOE) with regression modeling. By systematically processing experimental data, RSM elucidates the functional relationships between multiple independent input variables ( $X_i$ ) and the corresponding response variables (y). The primary objective of RSM is to characterize the interaction effects among various control factors. This facilitates the construction of an empirical model that predicts the response surface, thereby enabling the identification of the optimal process configuration for maximized performance and efficiency.

RSM of carbon emissions

Based on the preliminary ANOVA conducted in Section 3.2.1, this study identified the dominant processing parameters to streamline the optimization process. Parameters with negligible statistical influence—specifically bed temperature and print speed—were omitted from the subsequent RSM study to enhance model parsimony and predictive precision.

Consequently, three significant control factors were established: X₁ (layer height), X₂ (nozzle temperature), and X₃ (infill density). The coded and actual level settings for these variables are detailed in Table 10. To characterize the response behavior, a central composite design (CCD) comprising 20 experimental runs was implemented (Table 11). A second-order polynomial regression model for carbon emissions was developed using Minitab software to quantify the functional relationship between the input parameters and the environmental response. The estimated regression coefficients and corresponding statistical metrics are summarized in Table 12.

Table 10.

RSM level parameter design (carbon emissions).

Significant factor	Unit	-1	0	+1
$X_1$ : layer height	mm	0.15	0.2	0.25
$X_2$ : nozzle temperature	°C	190	200	210
$X_3$ : infill density	%	30	50	70

RSM: response surface method.

Table 11.

Design results of response surface method (carbon emissions).

Exp.	Layer height (mm)	Nozzle temperature (°C)	Infill density (%)	Carbon emissions (kg  ${\mathrm{CO}}_2$ e)
1	0.20	200.00	50.00	33.5
2	0.25	210.00	70.00	39.1
3	0.20	200.00	50.00	33.3
4	0.20	200.00	50.00	33.6
5	0.15	190.00	30.00	27.8
6	0.28	200.00	50.00	29.2
7	0.20	200.00	16.36	24.1
8	0.12	200.00	50.00	43.8
9	0.20	200.00	50.00	33.4
10	0.15	210.00	70.00	48.8
11	0.20	200.00	83.64	50.5
12	0.15	190.00	70.00	45.1
13	0.20	183.18	50.00	30.1
14	0.15	210.00	30.00	31.3
15	0.25	190.00	30.00	23.9
16	0.25	210.00	30.00	26.5
17	0.25	190.00	70.00	37
18	0.20	200.00	50.00	33.2
19	0.20	200.00	50.00	33.7
20	0.20	216.82	50.00	36.9

Table 12.

Regression model coefficients (carbon emissions).

Item	Effect	Coefficient	Coefficient standard deviation	T Value	p Value
Constant	–	33.481	0.364	92.02	0.000
X 1	−12.574	−6.287	0.406	−15.49	0.000
X 2	5.748	2.874	0.406	7.08	0.000
X 3	25.836	12.918	0.406	31.82	0.000
X 1 × X1	4.942	2.471	0.665	3.72	0.004
X 2 × X2	−1.058	−0.529	0.665	−0.8	0.445
X 3 × X3	6.542	3.271	0.665	4.92	0.001
X 1 × X2	−1.768	−0.884	0.892	−0.99	0.345
X 1 × X3	−6.434	−3.217	0.892	−3.61	0.005
X 2 × X3	−0.212	−0.106	0.892	−0.12	0.908

S = 0.892, R² = 99.27%, R² (adj.) = 98.61%, R² (pred.) = 94.53%.

Statistical significance was evaluated based on the p values summarized in Table 12, where values below 0.05 indicate a significant impact on the experimental response. The analysis indicated that the interaction effect between layer height ( $X_1$ ) and infill density (X₃) significantly influenced the carbon emissions response, evidenced by a p value of 0.005.

To quantify the functional relationship between the processing parameters and the environmental impact, an empirical second-order regression model was established using least squares. The resulting quadratic equation, expressed in terms of actual (physical) factors, was derived as follows:

Y = −96.8 − 33 $X_1$  + 1.053 $X_2$  + 0.360 $X_3$  + 349.5 $X_1\times X_1$  − 0.00187 $X_2\times X_2$  + 0.002891 $X_3\times X_3$ − 0.625 $X_1\times X_2$  − 1.137 $X_1\times X_3$  − 0.00019 $X_2\times X_3$

Where Y represents the estimated carbon emissions (kg ${\mathrm{CO}}_2$ e), $X_1$ denotes the layer height (mm), $X_2$ is the nozzle temperature (°C), and $X_3$ signifies the Infill density (%). The model demonstrates a high degree of fit, with an $R^2$ of 99.27%, indicating that the regression equation can account for the vast majority of the variance observed in the experimental data.

Figure 11 illustrates the 3D response surface plots generated via Minitab, depicting the functional relationships and interactions among the processing parameters. As evidenced by the plots, infill density ( $X_3$ ) was the most influential determinant of the environmental response, characterized by the steepest gradient along its axis. This high sensitivity indicates that carbon emissions increased substantially as the structural density of the specimen increased, primarily due to the prolonged extrusion time and increased material consumption.

Figure 11.

Factor contour diagram (carbon emissions): (a) surface plot of layer height and nozzle temperature, (b) surface plot of layer height and infill density, and (c) surface plot of nozzle temperature and infill density.

Conversely, layer height ( $X_1$ ) exhibited a significant inverse relationship with carbon emissions. Specifically, an increase in layer height reduced energy consumption, as a thicker layer decreased the total number of deposited layers and the associated travel path of the nozzle. In contrast, nozzle temperature ( $X_2$ ) displayed the least impact within the investigated experimental range, exhibiting only a subtle linear trend without significant curvature.

Furthermore, the plots confirmed a significant interaction effect between layer height and infill density ( $X_1$  ×  $X_3$ ), as evidenced by the distinct twisting and non-parallel nature of the surface profile. This visual observation corroborated the statistical significance (p = 0.005) reported in Table 12. Such an interaction implies that the influence of infill density on carbon emissions was modulated by the selected layer height, necessitating a holistic optimization approach rather than individual parameter adjustments.

In summary, to satisfy the “Smaller-the-Better” (STB) objective for minimizing carbon emissions, the optimal process configuration was located within the region characterized by low infill density, high layer height, and low nozzle temperature.

RSM of ultimate tensile strength

Based on the ANOVA results for UTS detailed in Section 3.2.2, bed temperature and print speed were confirmed as factors with lower influence following a two-stage statistical analysis. To streamline the optimization process, these parameters were excluded from the subsequent RSM study. The significant factors selected for refined optimization include $X_1$ (layer height), $X_2$ (nozzle temperature), and $X_3$ (infill density). Prior S/N ratio analysis indicated that the optimal solutions for $X_3$ and $X_2$ reside in the high-level range, while $X_1$ is located in the low-level range. Consequently, the RSM experimental levels were established to target the identified “potential optimal regions” at (0.10, 0.15, 0.20 mm), (200 °C, 210 °C, 220 °C), and (50%, 70%, 90%), respectively, as detailed in Table 13. The primary objective was to construct a predictive regression model within these regions to determine the absolute optimal parameter configuration. Following a central composite design (CCD) comprising 20 experimental runs (Table 14), a second-order polynomial regression model for UTS was developed utilizing Minitab. The resulting regression coefficients, quantifying the linear, quadratic, and interaction effects of the processing parameters, are summarized in Table 15.

Table 13.

RSM level parameter design (UTS).

Significant factor	Unit	−1	0	+1
$X_1$ : layer height	mm	0.10	0.15	0.20
$X_2$ : nozzle temperature	°C	200	210	220
$X_3$ : infill density	%	50	70	90

RSM: response surface method.

Table 14.

Design results of response surface method (UTS).

Exp.	Layer height (mm)	Nozzle temperature (°C)	Infill density (%)	Ultimate tensile strength (MPa)
1	0.20	200.00	90.00	23.97
2	0.15	210.00	70.00	30.89
3	0.20	220.00	90.00	25.45
4	0.10	200.00	50.00	34.5
5	0.15	210.00	70.00	30.89
6	0.15	210.00	36.36	26.16
7	0.20	220.00	50.00	21.21
8	0.23	210.00	70.00	17.71
9	0.07	210.00	70.00	44.06
10	0.10	220.00	50.00	36.63
11	0.15	193.18	70.00	29.33
12	0.15	210.00	70.00	30.89
13	0.15	226.82	70.00	32.44
14	0.15	210.00	70.00	30.89
15	0.10	220.00	90.00	43.96
16	0.15	210.00	70.00	30.89
17	0.15	210.00	70.00	30.89
18	0.20	200.00	50.00	19.97
19	0.10	200.00	90.00	41.41
20	0.15	210.00	103.64	35.1

Table 15.

Regression model coefficients (UTS).

Item	Effect	Coefficient	Coefficient standard deviation	T Value	p Value
Constant	(N/A)	30.887	0.104	297.48	0.000
X 1	−27.146	−13.573	0.116	−117.15	0.000
X 2	3.11	1.555	0.116	13.43	0.000
X 3	9.24	4.62	0.116	39.88	0.000
X 1 × X1	0.1	0.05	0.19	0.27	0.796
X 2 × X2	0.1	0.05	0.19	0.27	0.796
X 3 × X3	−0.41	−0.205	0.19	−1.08	0.306
X 1 × X2	−1.386	−0.693	0.255	−2.72	0.021
X 1 × X3	−4.242	−2.121	0.255	−8.33	0.000
X 2 × X3	0.466	0.233	0.255	0.92	0.381

S = 0.254579, R² = 99.94%, R² (adj.) = 99.88%, R² (pred.) = 99.51%.

To quantify the mechanical response of the 3D-printed specimens, the statistical significance of the processing parameters was evaluated using the p values summarized in Table 15. Factors with p values below 0.05 are considered to exert a statistically significant influence on the UTS. Notably, the interaction between layer height ( $X_1$ ) and infill density ( $X_3$ ) was identified as a critical determinant of the mechanical performance.

Using the method of least squares, a second-order polynomial regression model was established to describe the functional relationship between the input variables and the UTS. The resulting empirical equation, expressed in terms of actual factors, was formulated as follows:

Y = 16.0 − 8.1 $X_1$  + 0.062 $X_2$  + 0.1885 $X_3$  + 7.1 $X_1\times X_1$  + 0.000178 $X_2\times X_2$  − 0.000181 $X_3\times X_3$  − 0.490 $X_1\times X_2$  − 0.7500 $X_1\times X_3$  + 0.000413 $X_2\times X_3$

In this model, $Y$ represents the estimated UTS (MPa), while $X_1$ , $X_2$ , and $X_3$ denote layer height (mm), nozzle temperature (°C), and infill density (%), respectively.

The 3D response surface plots illustrated in Figure 12 show the interaction effects and parametric trends influencing mechanical performance. Analysis of these plots indicated that infill density ( $X_3$ ) and nozzle temperature ( $X_2$ ) were positively correlated with the UTS. Specifically, maximizing the infill density to 100% enhanced the internal structural robustness and increased the effective load-bearing area of the specimen. Simultaneously, elevated nozzle temperatures (up to 220 °C) reduced the melt viscosity of the polymer, thereby facilitating improved molecular diffusion and interlayer fusion across the deposition interfaces. Conversely, layer height ( $X_1$ ) exhibited a significant negative correlation with tensile strength. Adopting a minimum layer height of 0.10 mm consistently improved mechanical performance by refining the microstructure, increasing the contact area between adjacent strands, and minimizing internal porosity.

Figure 12.

Factor contour diagram (UTS): (a) surface plot of layer height and nozzle temperature, (b) surface plot of layer height and infill density, and (c) surface plot of nozzle temperature and infill density.

The significant interaction identified between layer height and infill density ( $X_1\times X_3$ ) further underscores the necessity of a holistic optimization approach, as the influence of one parameter was modulated by the level of the other. As observed in the twisting of the response surfaces, the optimal mechanical integrity was achieved through specific parameter combinations rather than independent adjustments. In summary, to attain the maximum UTS within the investigated experimental range, the optimal process configuration was characterized by a high infill density (100%), an elevated nozzle temperature (220 °C), and a refined layer height (0.10 mm). This configuration ensured a cohesive and robust internal architecture, thereby satisfying the “Larger-the-Better” objective for mechanical performance.

Gray relational analysis (GRA)

In the field of additive manufacturing (AM), reconciling reductions in carbon emissions with improvements in mechanical performance has emerged as a critical multi-objective optimization challenge. While the conventional Taguchi method is highly effective for optimizing individual parameters, its traditional application is inherently constrained by its focus on a single quality characteristic. In scenarios where performance metrics are mutually exclusive—such as the conflicting requirements between energy efficiency and structural integrity—the optimal factor levels for one objective often compromise the other. To resolve this trade-off and achieve a Pareto-optimal balance, this study integrated the Taguchi method with GRA to identify a globally optimal parameter combination that simultaneously satisfied multiple quality requirements.

The implementation of the GRA framework was executed through four systematic procedures:

• Step 1: data preprocessing: The raw experimental output data were normalized to a dimensionless range of 0–1 to eliminate the effects of differing units and scales.

• Step 2: ideal reference sequence selection: A standard reference sequence (typically set to 1) was established to represent the theoretical optimal target point for each response variable.

• Step 3: calculation of the gray relational coefficient (GRC): This coefficient quantifies the degree of correlation or proximity between the normalized experimental results and the ideal target sequence.

• Step 4: derivation of gray relational grade (GRG): The GRG was calculated as the weighted average of the GRC values, serving as a comprehensive multi-objective performance indicator.

A higher GRG value indicated that the corresponding parameter combination was closer to the ideal performance. This integrated Taguchi-GRA methodology has been applied effectively across diverse engineering applications, including the optimization of surface grinding parameters ⁴¹ and thermal management in immersion cooling systems. ⁴² The theoretical foundations and mathematical derivations employed in this analysis were based on established literature.^43,44 The specific computational results are detailed as follows.

Data pre-processing

The objective of data pre-processing is to normalize the raw experimental data for carbon emissions and UTS into a dimensionless range of [0, 1], thereby ensuring comparability across datasets with differing units and physical magnitudes. In this study, carbon emissions were treated as a “Smaller-the-Better” (STB) quality characteristic, given the primary environmental goal of minimizing the carbon footprint associated with the FDM process. To achieve this, the raw output was normalized using the following formula (equation (6)). Conversely, UTS was treated as a “Larger-the-Better” (LTB) characteristic to ensure superior mechanical performance of the printed components. This characteristic was processed using equation (7). The normalized experimental results for all 20 runs, derived from these transformations, are summarized in Table 16.

$$x_{\mathrm{ij}}^{\mathrm{carbon}}=\frac{max\left(y_{\mathrm{ij}}\right)-y_{\mathrm{ij}}}{max\left(y_{\mathrm{ij}}\right)-min\left(y_{\mathrm{ij}}\right)}$$

(6)

$$x_{\mathrm{ij}}^{\mathrm{UTS}}=\frac{y_{\mathrm{ij}}-min\left(y_{\mathrm{ij}}\right)}{max\left(y_{\mathrm{ij}}\right)-min\left(y_{\mathrm{ij}}\right)}$$

(7)

Where $y_{\mathrm{ij}}$ represents the original experimental sequence, $x_{\mathrm{ij}}$ denotes the sequence after data pre-processing, and $max\left(y_{\mathrm{ij}}\right)$ and $min\left(y_{\mathrm{ij}}\right)$ are the maximum and minimum values of the original sequence, respectively.

Table 16.

Data pre-processing.

Exp.	A	B	C	D	E	Carbon emissions ( $y_{\mathrm{ij}}$ )	$X_{\mathrm{carbon}}$	UTS ( $y_{\mathrm{ij}}$ )	$X_{\mathrm{UTS}}$
1	1	1	1	1	1	24.7	0.939	32.7	0.176
2	1	2	2	2	2	34.4	0.593	38.7	0.634
3	1	3	3	3	3	43.8	0.257	42.6	0.931
4	1	4	4	4	4	51	0	43.5	1
5	2	1	2	3	4	38.7	0.439	39.7	0.71
6	2	2	1	4	3	35.8	0.543	38.3	0.603
7	2	3	4	1	2	34.7	0.582	36.7	0.481
8	2	4	3	2	1	32.8	0.65	41.8	0.87
9	3	1	3	4	2	28	0.821	33.5	0.237
10	3	2	4	3	1	25.4	0.914	30.4	0
11	3	3	1	2	4	37.8	0.471	41.6	0.855
12	3	4	2	1	3	38.8	0.436	39.8	0.718
13	4	1	4	2	3	27.3	0.846	35.3	0.374
14	4	2	3	1	4	34.6	0.586	38.6	0.626
15	4	3	2	4	1	23	1	31.5	0.084
16	4	4	1	3	2	24.2	0.957	33.2	0.214

Determining the optimal target sequence

After the data normalization phase, the deviation sequence was computed to quantify the numerical difference between the processed experimental values and the ideal performance targets. In the GRA framework, the ideal reference sequence ( $x_0\left(k\right)$ ) was established as 1, representing the theoretical optimum for each quality characteristic. The deviation sequence, denoted as ${\mathrm{\Delta }}_{\mathrm{ij}}$ (k), was defined as the absolute magnitude of the difference between the reference sequence ( $x_0\left(k\right)$ ) and the normalized data sequence ( $x_{\mathrm{ij}}\left(k\right)$ ). This mathematical relationship was formulated in equation (8). The calculated deviation sequences for all 20 experimental runs are summarized in Table 17. These values were required for determining the global maximum deviation ( ${\mathrm{\Delta }}_{\max }$ ) and minimum deviation ( ${\mathrm{\Delta }}_{\min }$ ) within the entire dataset, which serve as foundational parameters for the subsequent calculation of the GRC.

$${\mathrm{\Delta }}_{\mathrm{ij}}\left(\mathrm{k}\right)=|x_0\left(k\right)-x_{\mathrm{ij}}\left(k\right)|=|1-x_{\mathrm{ij}}\left(k\right)|$$

(8)

Where:

  ${\mathrm{\Delta }}_{\mathrm{ij}}$ (k): The deviation sequence for the ith experiment and kth response,

  $x_0\left(k\right)$ : The ideal reference value (set to 1.0),

  $x_{\mathrm{ij}}\left(k\right)$ : The normalized experimental value derived from equation (7) or (8).

Table 17.

Deviation sequence results.

Exp.	A	B	C	D	E	$X_{\mathrm{carbon}}$	${\mathrm{\Delta }}_{\mathrm{carbon}}$	$X_{\mathrm{UTS}}$	${\mathrm{\Delta }}_{\mathrm{UTS}}$
1	1	1	1	1	1	0.939	0.061	0.176	0.824
2	1	2	2	2	2	0.593	0.407	0.634	0.366
3	1	3	3	3	3	0.257	0.743	0.931	0.069
4	1	4	4	4	4	0	1	1	0
5	2	1	2	3	4	0.439	0.561	0.71	0.29
6	2	2	1	4	3	0.543	0.457	0.603	0.397
7	2	3	4	1	2	0.582	0.418	0.481	0.519
8	2	4	3	2	1	0.65	0.35	0.87	0.13
9	3	1	3	4	2	0.821	0.179	0.237	0.763
10	3	2	4	3	1	0.914	0.086	0	1
11	3	3	1	2	4	0.471	0.529	0.855	0.145
12	3	4	2	1	3	0.436	0.564	0.718	0.282
13	4	1	4	2	3	0.846	0.154	0.374	0.626
14	4	2	3	1	4	0.586	0.414	0.626	0.374
15	4	3	2	4	1	1	0	0.084	0.916
16	4	4	1	3	2	0.957	0.043	0.214	0.786

The deviation sequences presented in Table 17 illustrate the performance extremes and the trade-offs in the experimental design. Specifically, Exp. 15 achieved a minimum deviation of zero for carbon emissions, indicating that this configuration reached the ideal performance benchmark for the environmental objective. Conversely, Exp. 4 yielded a minimum deviation of zero for UTS, identifying it as the optimal run for mechanical integrity.

These contrasting outcomes demonstrated the conflicting nature of the two optimization objectives: a parameter set optimized for sustainability often compromises structural performance, while a configuration optimized for maximum strength often results in a higher environmental impact. This tradeoff necessitated the subsequent calculation of the GRC to identify a robust trade-off point (or Pareto-optimal balance) that satisfied both mechanical and environmental requirements simultaneously.

Gray relational coefficient (GRC)

Upon determining the deviation sequence ${\mathrm{\Delta }}_{\mathrm{ij}}(k)$ , the GRC was computed to quantify the degree of correlation between individual experimental results and the ideal reference target. While the deviation sequence primarily represents an absolute numerical gap, the GRC transforms this distance into a relational closeness index within the interval [0, 1]. A GRC value of 1.0 signifies that the experimental result is perfectly consistent with the theoretical optimum, whereas a value closer to zero indicates a significant divergence.

The GRC ( ${\xi }_{\mathrm{iJ}}\left(k\right))$ was calculated according to equation (9), which incorporates a distinguishing coefficient to adjust the contrast between deviation values:

$${\xi }_{\mathrm{ij}}\left(k\right)=\frac{{\mathrm{\Delta }}_{\min }+\psi {\mathrm{\Delta }}_{\max }}{{\mathrm{\Delta }}_{\mathrm{ij}}(k)+\psi {\mathrm{\Delta }}_{\max }}$$

(9)

Where:

  ${\xi }_{\mathrm{iJ}}\left(k\right)$ : The gray relational coefficient for the jth experiment and kth response.

  $\psi$ : The distinguishing coefficient, typically set at 0.5 to provide moderate contrast and stability in the relational analysis.

  ${\mathrm{\Delta }}_{\min }$ and ${\mathrm{\Delta }}_{\max }$ : The global minimum and maximum deviations observed across all response characteristics in the entire experimental matrix.

In this formulation, $\psi$ represents the distinguishing coefficient, a pivotal parameter that modulates the mathematical resolution and contrast of the relational analysis. From a physical standpoint, this parameter attenuates the disproportionate influence of extreme deviation values, thereby enhancing the stability and robustness of the resulting GRCs. By tuning $\psi$ , the sensitivity of the GRC to inherent data fluctuations can be effectively regulated.

The value of the distinguishing coefficient is bounded within the interval $\psi \in$  [0, 1]. Adhering to established academic conventions in engineering optimization, a standard value of $\psi$  = 0.5 was adopted in this study. This specific value provides a moderate and unbiased distinguishing effect, ensuring an optimal balance between discriminative resolution and numerical stability without artificially skewing the results toward extreme deviations. Consequently, this theoretically grounded selection provided a rigorous foundation for multi-objective decision-making regarding carbon emissions and mechanical performance. The calculated coefficients for carbon emissions and UTS are summarized in Table 18.

Table 18.

GRC for carbon emissions and UTS.

Exp.	A	B	C	D	E	${\mathrm{\Delta }}_{\mathrm{carbon}}$	${\xi }_{\mathrm{carbon}}$	${\mathrm{\Delta }}_{\mathrm{UTS}}$	${\xi }_{\mathrm{UTS}}$
1	1	1	1	1	1	0.061	0.892	0.824	0.378
2	1	2	2	2	2	0.407	0.551	0.366	0.575
3	1	3	3	3	3	0.743	0.402	0.069	0.873
4	1	4	4	4	4	1	0.333	0	1
5	2	1	2	3	4	0.561	0.471	0.29	0.633
6	2	2	1	4	3	0.457	0.522	0.397	0.557
7	2	3	4	1	2	0.418	0.545	0.519	0.491
8	2	4	3	2	1	0.35	0.588	0.13	0.794
9	3	1	3	4	2	0.179	0.737	0.763	0.396
10	3	2	4	3	1	0.086	0.854	1	0.333
11	3	3	1	2	4	0.529	0.486	0.145	0.771
12	3	4	2	1	3	0.564	0.47	0.282	0.639
13	4	1	4	2	3	0.154	0.765	0.626	0.444
14	4	2	3	1	4	0.414	0.547	0.374	0.572
15	4	3	2	4	1	0	1	0.916	0.353
16	4	4	1	3	2	0.043	0.921	0.786	0.378

The results summarized in Table 18 provide a quantitative representation of the performance trade-offs inherent in the FDM process. Notably, Exp. 15 achieved a GRC of 1.000 for carbon emissions, signifying that this specific configuration minimizes the environmental footprint across the entire experimental design. However, the corresponding GRC for UTS was 0.353, indicating a substantial reduction in structural integrity.

Conversely, Exp. 4 represents the opposite extreme, attaining a GRC of 1.000 for UTS but yielding a GRC for carbon emissions of only 0.333. This confirms that the parameters required for maximum mechanical strength—typically high infill density and elevated temperatures—led to the highest environmental impact within the tested range.

These extreme cases demonstrated the inherent limitations of single-objective optimization, where pursuing the theoretical optimum for one quality characteristic inevitably degrades the other. To resolve this multi-objective dilemma, the gray relational grade (GRG) was calculated as a weighted average of these coefficients. This integrated metric facilitated the identification of a global optimal solution (or Pareto-optimal point) that effectively reconciles the requirements for both environmental sustainability and mechanical performance.

Gray relational grade (GRG)

The GRG, denoted as ${\gamma }_i$ , serves as the definitive metric for the multi-objective optimization phase. Its primary function is to consolidate the discrete GRC values of multiple individual objectives into a singular, comprehensive index. This procedure effectively transforms a multi-dimensional optimization challenge into a single-objective evaluation, quantifying the integrated performance of each experimental run across all conflicting indicators.

The GRG is derived by calculating the arithmetic mean of the GRC values for each experimental group, as formulated in equation (10):

$${\gamma }_i=\frac{1}{m}\sum _{k=1}^m{\xi }_i\left(k\right)$$

(10)

Where:

  ${\gamma }_i$ : The gray relational grade for the ith experimental run.

  ${\xi }_i(k)$ : The gray relational coefficient for the kth response characteristic in the ith run.

 m: The total number of performance indicators (in this study, m = 2, representing UTS and carbon emissions).

A higher GRG value indicates that the collective performance of the respective parameter combination exhibits superior proximity to the ideal reference sequence. Within the context of this study, a maximum GRG represents a balanced performance between mechanical strength and environmental sustainability. Consequently, the parameter combination yielding the highest GRG is designated as the global optimal solution for the multi-objective problem. The comprehensive GRG values and their corresponding rankings for all 20 experimental runs are summarized in Table 19.

Table 19.

GRG and ranking for each experimental run.

Exp.	A	B	C	D	E	${\xi }_{\mathrm{carbon}}$	${\xi }_{\mathrm{UTS}}$	GRG	Rank
1	1	1	1	1	1	0.892	0.378	0.635	6
2	1	2	2	2	2	0.551	0.575	0.564	11
3	1	3	3	3	3	0.402	0.873	0.641	5
4	1	4	4	4	4	0.333	1	0.667	3
5	2	1	2	3	4	0.471	0.633	0.552	14
6	2	2	1	4	3	0.522	0.557	0.54	15
7	2	3	4	1	2	0.545	0.491	0.518	16
8	2	4	3	2	1	0.588	0.794	0.691	1
9	3	1	3	4	2	0.737	0.396	0.566	10
10	3	2	4	3	1	0.854	0.333	0.593	9
11	3	3	1	2	4	0.486	0.771	0.631	7
12	3	4	2	1	3	0.47	0.639	0.554	13
13	4	1	4	2	3	0.765	0.444	0.605	8
14	4	2	3	1	4	0.547	0.572	0.559	12
15	4	3	2	4	1	1	0.353	0.677	2
16	4	4	1	3	2	0.921	0.378	0.655	4

As demonstrated by the results in Table 19, Exp. 8 attained the peak GRG, thereby being identified as the global optimal parameter combination within the investigated experimental matrix. The superior performance of this run indicated its proximity to the ideal reference sequence, effectively reconciling the divergent requirements of environmental sustainability and structural integrity.

The optimal parametric configuration was identified as A2B4C3D2E1, which corresponds to the following processing levels: layer height of 0.15 mm, nozzle temperature of 220 °C, bed temperature of 70 °C, print speed of 50 mm/s, and infill density of 30%.

This combination represents a balanced configuration for FDM processing, as it maintains relatively low carbon emissions without a substantial reduction in UTS. From an industrial engineering perspective, these findings provide highly valuable reference data for manufacturers seeking to implement sustainable manufacturing protocols while maintaining high-performance standards. Adoption of this optimized configuration resulted in a significant reduction in the environmental footprint while preserving the mechanical functionality of the printed components.

Verification analysis of the optimal parameter combination

Performance comparison: Multi-objective versus single-objective optimization

This section evaluates the efficacy of the GRA in reconciling the inherent conflict between UTS and the carbon emissions in the FDM process. To validate the robustness of the multi-objective optimal configuration (A2B4C3D2E1), its performance was compared against two boundary scenarios derived from single-objective optimization:

This section evaluates the efficacy of the GRA in reconciling the inherent conflict between UTS and the carbon emissions in the FDM process. To validate the robustness of the multi-objective optimal configuration (A2BC3D2E1), its performance was compared against two boundary scenarios derived from single-objective optimization:

1. Minimum carbon emission configuration (A4B1C1D3E1).

2. The maximum UTS configuration (A1B4C3D2E4).

The detailed performance indicators and the calculated Improvement Rates (IR) are summarized in Tables 20 and 21. To provide a quantitative assessment of the optimization gains, the IR for carbon emissions was calculated using the “Smaller-the-Better” (STB) characteristic (equation (11)), while the IR for UTS was determined using the “Larger-the-Better” (LTB) characteristic (equation (12)):

$${\mathrm{IR}}_{\mathrm{STB}}=\frac{(X_{\mathit{Single}}-X_{\mathrm{dual}})}{X_{\mathit{Single}}}\times 100\%$$

(11)

$${\mathrm{IR}}_{\mathrm{LTB}}=\frac{(X_{\mathrm{dual}}-X_{\mathit{Single}})}{X_{\mathit{Single}}}\times 100\%$$

(12)

Where $X_{\mathit{Single}}$ and $X_{\mathrm{dual}}$ represent the results from single-objective and multi-objective (dual-target) optimized configurations, respectively.

Table 20.

Comparison between multi-objective optimization and carbon emission-oriented optimization.

	Optimal parameter combination for dual objectives	Optimal parameter combination for carbon emissions	Percentage improvement
Optimal combination	A2B2C3D2E1	A4B1C1D3E1
Carbon emissions	35.5	22.9	-55.02%
Ultimate tensile strength (MPa)	38.1	32.1	+18.69%

Table 21.

Comparison between multi-objective optimization and UTS-oriented optimization.

	Optimal parameter combination for dual objectives	Optimal parameter combination for ultimate tensile strength	Percentage improvement
Optimal combination	A2B2C3D2E1	A1B4C3D2E4
Carbon emissions	35.5	42.2	15.88%
Ultimate tensile strength (MPa)	38.1	44.2	-13.80%

Based on the empirical data in Tables 20 and 21, prioritizing a single quality characteristic resulted in performance trade-offs at the extremes. The configuration optimized solely for minimum carbon emissions (A4B1C1D3E1) significantly compromises mechanical integrity, yielding a UTS of only 32.1 MPa. Conversely, maximizing UTS (A1B4C3D2E4) imposed the maximum environmental burden within the tested range, with carbon emissions reaching 42.2 kg ${\mathrm{CO}}_2$ .

In contrast, the integrated GRA approach identified a synergistic balance between these mutually exclusive targets. The derived optimal combination (A2B2C3D2E1) effectively resides on the Pareto Front, successfully mitigating the sharp performance drops observed in the boundary configurations.

As quantified by the improvement metrics:

• Mechanical integrity: The balanced solution provided a UTS, that is, 18.69% above the minimum carbon emission set.

• Environmental sustainability: Simultaneously, it achieved a 15.88% reduction in carbon emissions compared to the maximum UTS set.

These results indicate that the GRA-based configuration reduced the “fragility” of single-objective models. By sacrificing a limited portion of the theoretical maximum in one dimension, the system achieved improved viability in the other, thereby providing a more balanced solution for engineering applications requiring both mechanical performance and environmental considerations.

Comparative analysis of optimization methodologies

To evaluate the efficacy of different optimization frameworks, a benchmark comparison between GRA, the Taguchi method, and RSM is presented in Table 22. This comparison illustrates the performance differences between multi-objective and single-objective optimization strategies.

• Carbon emission optimization: As shown in Table 22, the absolute minimum carbon emission of 22.9 was obtained using the single-objective parameter combination A4B1C1D3E1, identified by both the Taguchi method and RSM. While the dual-objective GRA did not reach this absolute minimum, it produced the second-lowest emission value, maintaining a competitive environmental profile.

• Mechanical performance (UTS): Regarding structural integrity, the peak UTS was recorded at 44.2 MPa. This was attained through the RSM single-objective optimal combination A1B4C3D2E4. In alignment with the environmental trends, the dual-objective GRA results ranked second in this category.

Figure 13 illustrates the trade-offs and integrated performance of the experimental models across three scenarios: (a) minimum carbon emission combination; (b) maximum UTS combination; and (c) dual-objective GRA optimal combination.

Table 22.

Comparison of the GRA and the Taguchi method.

Control factor	Optimal parameter combination for carbon emissions	Optimal parameter combination for ultimate tensile strength (MPa)	Optimal parameter combination for dual objectives
A. Layer height (mm)	0.25	0.10	0.15
B. Nozzle temperature (°C)	190	220	200
C. Bed temperature (°C)	50	70	70
D. Print speed (mm/s)	60	50	50
E. Infill density (%)	30	90	30
Carbon emissions	22.9	42.2	35.5
Ultimate tensile strength (MPa)	32.1	44.2	38.1

Figure 13.

Physical specimens and fracture morphology of 3D-printed parts under different optimization strategies: (a) minimum carbon emission combination, (b) maximum UTS combination, and (c) dual-objective GRA optimal combination.

Conclusions

This study developed an optimization framework by integrating the Taguchi method, RSM, and GRA to address multi-objective challenges in the FDM 3D printing process. By reconciling the trade-off between environmental sustainability (minimized carbon footprint) and mechanical integrity (maximized UTS), this research offers a structured analytical basis for sustainable production strategies. The following conclusions were drawn:

1. Identification of predominant influencing factors: ANOVA results indicated that infill density and layer height are the most critical determinants of FDM process quality:

2. Ultimate tensile strength (UTS): Infill density remained the primary contributing factor (45.51%), followed by layer height (23.69%) and nozzle temperature (14.82%). Enhanced mechanical performance was associated with higher structural density and improved interlayer fusion facilitated by thinner layers and elevated temperatures.

Constraints of single-objective optimization: The investigation reveals that prioritizing a single quality characteristic resulted in detrimental performance extremes that are often non-viable in engineering practice:

The “environmental priority” configuration (A4B1C1D3E1) successfully minimized emissions to 22.9 but yielded an insufficient UTS of 32.1 MPa, compromising structural reliability.

The “mechanical priority” configuration (A1B4C3D2E4) achieved a peak UTS of 44.2 MPa but produced carbon emissions of 42.2, increasing environmental impact within the tested range.

Efficacy of the multi-objective GRA framework: The GRA effectively addressed the identified trade-off by consolidating conflicting responses into a single GRG. The global optimal parameter combination was determined as A2B2C3D2E1, characterized by the following settings: layer height of 0.15 mm, nozzle temperature of 200 °C, bed temperature of 70 °C, print speed of 50 mm/s, and infill density of 30%.

Verification and empirical gains: Verification experiments confirmed that the GRA-optimized configuration approximated a Pareto-efficient solution within the investigated parameter range:

Versus the minimum carbon benchmark: The balanced GRA solution increased UTS by 18.69% while also leading to a moderate increase in carbon emissions.

Versus maximum UTS benchmark: The GRA solution achieved a 15.88% reduction in carbon emissions while maintaining the mechanical strength within a high-performance threshold.

This research demonstrates that the Taguchi-GRA methodology functions as a robust Pareto-optimization framework for the industrial sector. While the absolute numerical findings are grounded in the specific thermodynamic profile of PLA and the kinetic constraints of the DUAL300 system, the underlying physical mechanisms—such as layer height acting as an energy-phase emission multiplier and infill density dictating material-phase emissions—represent universal characteristics inherent to all extrusion-based additive manufacturing processes.

Crucially, this framework acknowledges the inherent antagonism between mechanical integrity and environmental impact, providing a systematic, material-agnostic approach to navigating this controlled trade-off. By adopting the Pareto-optimal configuration identified in this study, manufacturers can achieve a verifiable 15.88% reduction in their carbon footprint while maintaining a safety factor that exceeds common structural requirements for functional industrial components. This optimization aligns with the practical reality that most 3D-printed parts do not operate at their theoretical breaking point, thereby offering a pragmatic pathway for the industry to reduce energy expenditures, align with global carbon neutrality mandates such as ISO 14064, ⁴⁵ and enhance Environmental, Social, and Governance (ESG) performance in diverse smart factory environments across various hardware and material platforms.