Lean-and-Green Datacentric Engineering in Laser Cutting: Non-Linear Orthogonal Multivariate Screening Using Gibbs Sampling and Pareto Frontier

Abstract

Metal processing may benefit from innovative lean-and-green datacentric engineering techniques. Broad process improvement opportunities in the efficient usage of materials and energy are anticipated (United Nations Sustainable Development Goals #9, 12). A CO₂ laser cutting method is investigated in this study in terms of product characteristics (surface roughness (SR)) and process characteristics (energy (EC) and gas consumption (GC) as well as cutting time (CT)). The examined laser cutter controlling factors were as follows: (1) the laser power (LP), (2) the cutting speed (CS), (3) the gas pressure (GP) and, (4) the laser focus length (F). The selected 10mm-thick carbon steel (EN10025 St37-2) workpiece was arranged to have various geometric configurations so as to simulate a variety of real industrial milling demands. Non-linear saturated screening/optimization trials were planned using the Taguchi-type L₉(3⁴) orthogonal array. The resulting multivariate dataset was treated using a combination of the Gibbs sampler and the Pareto frontier method in order to approximate the strength of the studied effects and to find a solution that comprises the minimization of all the tested process/product characteristics. The Pareto frontier optimal solution was (EC, GC, CT, SR) = (4.67 kWh, 20.35 Nm³, 21 s, 5.992 μm) for the synchronous screening/optimization of the four characteristics. The respective factorial settings were optimally adjusted at the four inputs (LP, CS, GP, F) located at (4 kW, 1.9 mm/min, 0.75 bar, +2.25 mm). The linear regression analysis was aided by the Gibbs sampler and promoted the laser power and the cutting speed on energy consumption to be stronger effects. Similarly, a strong effect was identified of the cutting speed and the gas pressure on gas consumption as well as a reciprocal effect of the cutting speed on the cutting time. Further industrial explorations may involve more intricate workpiece geometries, burr formation phenomena, and process economics.

1. Introduction

Laser cutting is an advanced manufacturing process that has found many applications in fabricating various types of products [1]. From drilling holes to milling sheet plates, laser cutting is recommended for superior surface integrity results. From a machining perspective, the process uniqueness of laser beam cutting is identified to deliver focused thermal energy to the targeted workpiece, without any machine tool-end coming in direct contact with the processed item. An additional advantage stems from the fact that this machining approach may even be used for treating hard steel alloys, which are known for their remarkable mechanical properties [2]. Specifically, modelling a CO₂ laser cutting process is a subject that still attracts a lot of attention. There is a plethora of mechanical, metallurgical, and geometrical aspects to be probed due to the growing number of innovative materials that permeate the diverse areas of modern product development. Accordingly, datacentric engineering techniques are deployed to facilitate the uncovering of improvement opportunities in highly customizable products, where often there is no prior experience. Gathering insights for fruitful problem-solving is benefited by an inclusive perspective that engages all available areas of theoretical manipulation, i.e., analytical, statistical, and algorithmical toolboxes. What makes the contemporary laser cutting processes more intriguing is the international initiative for the urgent greening of new products and technologies [3]. It is imperative for pollution curbing to be engrained in innovative production technologies in order to advance products that will remain sustainably useful in the future [4,5]. Resilient process design is greatly advantageous and certainly attainable by combining green and lean principles to lead successful endeavors in sustainability enhancement [6,7,8]. The lean aspect is an integral part of achieving better operational performance by minimizing production costs and waste creation [9,10,11,12]. Of course, lean-and-green practices play a central role in supporting the neoteric circular economy; it provides the impetus for zero waste generation and maximum exploitation of the resources from ‘cradle-to-grave’ [13,14,15,16,17]. To improve the environmental performance of complex manufacturing processes, new methodologies are necessitated to incorporate decision-making tools and techniques during the improvement project implementation. Discovering opportunities for greening the future digital factory requires high-level datacentric engineering techniques that seamlessly incorporate lean management techniques to the extensive quality improvement toolbox (Six Sigma) [18,19,20,21,22,23,24].

A lean-and-green study of a laser cutting process involves a collection of product and process characteristics. Product characteristics may involve any of the physical material characteristics such as surface roughness, surface dross, top and bottom kerf width, kerf taper, heat-affected-zone width, striation formation, and ablation depth. Lean-and-green process characteristics may be considered to be the machining-center energy consumption, the assist gas consumption, and the cutting time. The cutting time is a lean process characteristic because it explicitly represents the key productivity level of the laser cutting process which is directly associated to the process cycle time.

To invoke a process improvement study, it is necessary to consider the implicated laser beam factors such as the beam wavelength, the laser pulse power, the beam polarization, the beam type, and the power intensity. Laser machining process factors may be the standoff distance, the nozzle diameter, the gas type, the gas pressure, the laser beam focus, and the cutting speed. Furthermore, material types with different combinations of geometrical shapes and dimensions conclude the potential group of product characteristics. Based on the current research literature, there is no such comprehensive study that includes all the above-mentioned controlling factors. Similarly, no study has taken into account all the material process characteristics that were listed above. There is a reasonable explanation for this absence of all-inclusive investigations; it is related to the practicality of carrying out such studies, as well as to the associated costs of completing such a research effort. Despite the speed of cutting being high and the machining cost being low, the inherent high accuracy of laser beam cutting and its high level of production applications cannot compensate for the bulk of experiments that need to be programmed in order to carry out a comprehensive product/process improvement project. A substantial amount of time needs to be dedicated for the scheduled trials which should be conducted on the shop floor. Usually, this time must be subtracted from the regular production time of the machining center.

Along those lines, there is generally an increasing interest in sustainable machining that is eco-friendly, waste-free, and energy efficient [25]. Modeling the enhancement of the laser cutting performance, through a combination of screening and optimization, requires the deployment of advanced statistical methods, which should be selected to be flexible enough to steer the optimal solution according to the appropriately formulated and weighted product/process characteristic preferences [26,27]. It is clearly anticipated that realistic performance improving studies are constrained to deal with multi-objective screening/optimization schemes. For materials with innate superior performance, such as stainless steel, the multi-response approach appears to be the Grey Relational Analysis (GRA) [28,29]. Nevertheless, GRA is a rather subjective method because of the requirement of casting the problem in a single “master” response, often resorting to uniform coefficients of identification, which are introduced to weigh the implied responses. The resulting grey relational grade values are converted to ranks. A weakness of the method is that no multifactorial non-parametrics are properly applied to gauge the differences among the factorial settings. An alternative method is to keep the multi-variate dataset small, using the Taguchi methods, and then combine neural networks and data envelopment analysis to obtain a desirable solution [30]. Computer-aided Taguchi samplers are shown to require fuzzification of the multi-objective problem, before approximating a final prediction [31]. Taguchi methods may be assisted by the implementation of a response surface methodology in order to refine the optimal solution, in addition to providing a comparative measure on the proposed solution, in the case of investigating single characteristics [32]. Nevertheless, for a laser cutting study to be meaningful to the shop floor personnel, the multifactorial screening/optimization process cannot be restricted to single characteristics. For example, recent reviews on laser cutting techniques have shown that process-specific characteristics, such as the assist gas, cannot be excluded from a rigorous examination of a laser cutting process [33,34]. It is agreeable that robust treatments should be favored in the comprehensive assessments of laser cutting conditions and performances [35]. The laser cutting performance-improvement attempt should incorporate information from the surface integrity metrics that capture the metal workpiece micro-geometry in a convenient fashion which is also easily interpretable [36]. Utilizing more traditional cutting techniques, it was shown that lean-and-green indicators may be flexibly inserted in the broader multi-response screening/optimization scheme [37]. It was exemplified that weighted synchronous multifactorial-setting adjustments may be attained in a material removal process where the machined specimens consisted of an austenitic stainless steel product.

The novelty of this work is based on examining the performance of a laser cutting process by monitoring a key product indicator of surface integrity, the workpiece surface roughness. The showcased material is a carbon steel product, where different geometrical shapes are drawn on the designed item. However, the emphasis of this multivariate screening/optimization exercise is placed on the lean-and-green aspect of the laser cutting procedure. Three process characteristics will be tracked down, alongside the surface roughness: the energy consumption, the gas consumption, and the cutting time. Moreover, four controlling factors will be profiled for their relative influence on the laser cutting process quality. Their optimal adjustments take into account any potential non-linear trends. The examined controlling factors are the laser power, the gas pressure, the cutting speed, and the laser-beam focus length.

To economize the requisite trial volume, a non-linear family of fractional factorial designs (FFDs) is chosen that demands a minimum number of preset runs [38]. FFDs significantly condense the number of experiments, thus minimizing the research logistics requirements, besides their associated incurring costs. This convenience permits the utilization of production machinery to be reserved for testing an experimental product, in pragmatic conditions, in addition to minimizing production schedule interruptions owing to the committed processor. A desired trait of the FFD-generated data is that they may be straightforwardly analyzed by ordinary multifactorial statistical tools, such as analysis of variance (ANOVA) and General Linear Modelling (GLM) [39]. The particular FFD data-sampler is selected among the available non-linear Taguchi-type orthogonal arrays (OAs), because of its uniqueness in balancing its non-linear trial-run structure among all arranged factors [40,41]. The most economic trial planner that accommodates non-linear relationships into a compact design and serves as many as four controlling factors is the L₉(3⁴) OA. The replicated trial plan will be treated with two different approaches: (1) to profile the strong effects for each process/product characteristic, and (2) to obtain an optimal multifactorial solution by synchronously minimizing all four examined characteristics. To accomplish the former task, the Gibbs sampling method is used to form the multi-distribution marginal densities [42,43,44,45]. Consequently, a robust application of linear regression can be applied to size the factorial effects through the approximated coefficients of regression. To estimate an optimal solution, the Pareto frontier method is implemented [46,47,48,49,50,51,52,53,54]. The Pareto frontier approach has been used regularly in quality improvement studies that seek to simultaneously optimize several product and process characteristic performances [55,56,57,58]. The implemented Skyline queries and Hasse diagrams are especially practical means of prioritizing the multi-response data representation. They allow assessment of a solution by forming a distinct border line where the optimal solution is bound to be located.

The paper has been structured as follows: (1) a methodology section provides the technical details about the laser cutting procedure and the carbon steel material that is used in the experiments, as well as a quick theoretical overview of the employed statistical techniques, the Gibbs-sampling linear regression and the multi-objective Pareto frontier, along with a synopsis of the data science modules which are implemented on the R-platform, (2) in the Results section the datacentric profiling outcomes highlight the Gibbs-sampled coefficients of regression estimations, and the skyline and Hasse diagrams, furnishing suggested scenarios for the optimal multi-objective solutions, (3) a Discussion section comments on the statistical assumptions that were necessary to simplify some parts of the screening/optimization methodology, and (4) the key findings and future work are summed up in a Conclusion section.

2. Materials and Methods

2.1. The Experimental Setup Requirements

2.1.1. Description of the Material’s Specifications and Geometries

The selected material is a widely used carbon-and-steel alloy—EN10025 St37-2 plate (DIN17100). The investigated plate thickness was 10mm, which is characterized by a tensile strength (Rm transv.) of 340–470 N/mm², and a fracture elongation (transv. min) of 24%. The chemical element upper limits by mass are as follows: C = 0.17% (max), Mn = 1.40% (max), P = 0.045% (max), and S = 0.045% (max). The dimensioning drawing specifics of a part are shown in Figure 1. The particular geometry and cutting paths were determined such as to examine synchronously the cutting quality of the specimen by permitting the laser cutting of two differently shaped workpiece corners: a round corner and a sharp corner. To complete the cutting process challenge, a disk was to be removed from the center of the test plate.

2.1.2. Description of the Laser Cutting System

In this work, the CO₂ laser cutting machinery is a Trumpf TruLaser 5040 (Trumpf, Ditzingen, Germany) with maximum laser power of 6 kW, which is suitable for cutting sheets from 0.5 to 30 mm. The CO₂ laser wavelength is 10.600 nm (far-infrared). The auxiliary gas mixture consists of the elements N₂ and He. The mirror-deflected beam is focused with either a lens or a mirror. During laser cutting, the assist gas is blown out through a nozzle with its diameter for this group of experiments set at 1.2 mm. The size of the nozzle diameter is important in order to calculate the gas consumption in Nm³/h units, given that what is known in real time during the cutting process is the applied gas pressure in bars.

2.1.3. Assignment of Controlling Factors and Settings for the Laser Process

The examined laser process controls were as follows: (1) the laser power, (2) the cutting speed, (3) the gas pressure, and (4) the focus length. In

Table 1. The controlling factor list and their studied settings for the laser cutting experiments.

				Settings
Controlling Factors	Abbreviation	Units	Level 1	Level 2	Level 3
Laser Power	LP	W	4000	5000	6000
Cutting Speed	CS	m/min	0.8	1.6	2.4
Gas Pressure	GP	bar	0.5	0.7	1
Focus Length	F	mm	−4	+2.2	+6.5

, we list the four selected controlling factors, their abbreviations, their units, and their three tested adjustments. The experimental design was kept lean by implementing an L₉(3⁴) OA scheme [40,41] that automatically reduced the total volume of trials to just 1/9 of what would be expected from a full-factorial plan. Realizing that the complex phenomena behind the laser cutting process would require a non-linear modeling, the nine-run three-level OA was proposed to capture with a minimum effort (another lean opportunity) any underlying curvature effects. In

Table 2. The non-linear nine-run saturated L₉(3⁴) OA with their respective factorial assignments and their physical adjustments.

Trial No.	LP	CS	GP	F
1	4000	0.8	0.5	−4
2	4000	1.9	0.75	2.25
3	4000	3	1	6.5
4	5000	0.8	0.75	6.5
5	5000	1.9	1	−4
6	5000	3	0.5	2.25
7	6000	0.8	1	2.25
8	6000	1.9	0.5	6.5
9	6000	3	0.75	−4

, the L₉(3⁴) OA columns are assigned to each of the four controlling factors, and the matrix elements are filled with their physical settings from Table 1.

2.1.4. The Process/Product Characteristics in the Lean-and-Green Design

There are several laser performance characteristics that may provide all-around surface integrity of the processed part. In this work, only the surface roughness was measured, since it is an effective indicator of the quality status of the machined area of the workpiece. Surface roughness (SR) is an important performance characteristic because it allows the synchronous assessment of both, product and process, capabilities. Usually, the surface roughness is a product characteristic which is sought to be minimized. If the unit surface is assessed to be smooth for a particular application, then the achieved performance may be adequate as long as the measurements are shown to be repeatable. Generally speaking, for this type of processing, the surface is considered smooth with the proviso that the laser cutting performance is maintained at Ra < 6.3 μm. Surface roughness measurements (Figure 2) were carried out by a portable roughness gauge (CARLSURF 9J008, RugoSURF, TESA, Renens, Switzerland).

From a productivity perspective, the cutting time (CT) is a lean process characteristic that needs to be minimized. It is important because it directly affects unit cost and it determines the production capacity. The Trumpf TruLaser 5040 numerical controller allows programming of the machining task based on the inputs of the four controlling factors (Figure 3A). Therefore, the cutting time was measured in seconds, and it was recorded for the total batch of the scheduled workpieces (Figure 3B). The energy consumption (EC) is a lean-and-green process characteristic. It is a quantity that is sought to be minimized, and it is measured in kWh. It is a highly critical process characteristic because it significantly influences the unit product cost. The energy consumption was calculated from the power consumption of the laser cutting center and the cutting time; both are easily obtained. In the stand-by mode, the power for the laser cutting center was 23 kW. The laser power was preset according to the trial schedule of Table 2, and the relationship between the laser output as percentage of the maximum output and the power consumption of the machine and laser is available in the form of a graph [59]. The auxiliary gas consumption is also a process characteristic that was sought to be minimized, and it is also viewed as a lean-and-green indicator. It was calculated in a straightforward manner from the linear dependence of gas consumption (Nm³/h) on the applied gas pressure (bar), given that the slope of the line is regulated by the selected nozzle diameter. The calculation details may be found in ref. [59]. The recorded total gas consumption was converted in units of Nm³, after taking into account the duration for completing a processed specimen.

The trial replication tactic was also decided based on lean-and-green criteria. The experiments were initially scheduled for at least duplicating the OA trial runs. If the variation between two replicated trials were judged to be substantial, then additional rounds of replicates would have been conducted, until the overall variability of the replicated measurements was moderated. This way, it was thought that additional trial cost savings could have been attained by controlling the extent of replication. The lean-and-green aspect was recognized in the opportunity to further minimize the usage of resources, manhours, and equipment unavailability.

2.2. The Theoretical Modeling of the Multivariate Screening Optimization

The non-linear OA sampler which is proposed in this application is an FFD-based planner suitable for handling fixed factor effects in lean-and-green datacentric engineering applications. The OA-sampler designates the combination of the factorial levels which comprise a prescribed formulation matrix; each row is a unique factorial recipe. The formulation matrix is an m × n array, in which the n rows represent the experimental run specifics, and the m columns represent the considered controlling factors. The three-level OA schemes unite fixed-factor multi-parameter linear modeling with the possibility to uncover quadratic contributions. In saturation mode (maximum utilization), the three-level OAs render an expression between the number of runs and the number of trials: n = 2 × m + 1. The manipulated m controlling factors are defined as follows: X_j for 1 ≤ j ≤ m (m ϵ N), and their respective factor settings are denoted as x_ij for 1 ≤ i ≤ n (n ϵ N), and 1 ≤ j ≤ m. Non-linear OA schemes are constructed with k_j fixed levels for each jth factor (1 ≤ j ≤ m) and 3 ≤ k_j ≤ K_j (K_j ϵ N). The multi-characteristic response matrix R = {r_icd}, with 1≤ i ≤ n, 1 ≤ c ≤ L (L ϵ N), and 1 ≤ d ≤ D (D ϵ N) is generated by L characteristic responses, R_c; each c^th matrix column is replicated D total times.

2.2.1. The Pareto Frontier Profiler

The Pareto frontier solution is a set of Pareto efficient solutions for multi-objective screening optimization [46,47,48,49,50,51,52,53,54,55]. In a generalized manner, the multi-objective optimization problem is cast to the form minimize_X R(X) = [R₁(X), R₂(X), …, R_c(X)]^T according to conditions h_i(X) ≤ 0, with i = 1, 2, …, s being the number of equality and inequality constraints. A point x_ij* is Pareto optimal or weakly Pareto optimal if there does not exist another point x_ij such that R(X) ≤ R(X*) (and at least for one characteristic R_c(X) < R_c(X*)) or R(X) < R(X*), respectively. Thus, it is not certain that a single solution will satisfy all optimization goals for all criteria: min {R₁, R₂, R₃, …, R_L).

Even though a nondominated Pareto optimal solution exists, it may be contested that it is within a confidence interval optimal in case the R matrix elements are random variables. Therefore, the Pareto frontier is bounded by the Nadir objective vector: v_n = {sup_X*j R_c} ∀ 1 ≤ j ≤ m and 1 ≤ c ≤ L, and the ideal objective vector v_i = {inf_X*j R_c} ∀ 1 ≤ j ≤ m and 1 ≤ c ≤ L.

2.2.2. The Gibbs Sampling in Simple Regression Modeling

The basic general linear modeling of the multivariate characteristics R on the factorial settings X is simply given for the coefficients of regression vector β by:

R = Xβ + ε with ε ~ N(0, σ²) and β ~ N(μ, σ²) and σ²~Γ⁻¹(α,β)

Therefore, the posterior conditional probabilities for the coefficients of regression are conveniently expressed as:

p(β, σ²|R) ∝ N(Xβ, σ²) N(μ, τ²) Γ⁻¹(α,β)

and hence, the distribution of the coefficients of regression may be approximated as:

β|R, σ² ~ N(m, s) with m = (X′X)⁻¹ X′ R and s = σ² (X′X)⁻¹

where σ²|R, β ~ N(Xβ, σ²) Γ⁻¹(α,β) ~ Γ⁻¹(α_n,β_n)

given that α_n = a+ n/2 and β_n = (R − Xβ)′ (R − Xβ)/2

2.3. The Methodological Outline

The proposed methodology may be recapitulated as follows:

(1)

Determine the relevant cutting process characteristics that include product characteristics, i.e., surface roughness, but also lean-and-green process characteristics such as the energy and gas consumptions and cutting time.

(2)

Select a group of controlling factors pertinent to the laser beam cutting process.

(3)

Determine the operating end points for each controlling factor, from step 2, and ensure that there is adequate representation of factor settings to retrieve potential curvature effects.

(4)

Select an OA-sampler that best accommodates the factorial settings that were decided in step 3.

(5)

Carry out the minimum number of replications for the assigned fractional factorial recipes that were formulated from step 4.

(6)

Test the sufficiency of the extent of replication using linear regression methods for slope and intercept drifts.

(7)

Use boxplot-based response graphs to identify potential strong regressors.

(8)

Test process/product characteristics for two-way correlations, and decide to retain only those characteristics that provide unique information to the screening/optimization problem.

(9)

Combine the visual and numerical toolset for the Pareto frontier analysis to find the optimal factorial recipe by utilizing: (1) the Hasse diagrams and (2) the Pareto skylines.

(10)

Implement Gibbs sampling to obtain the posterior conditional distributions of the multivariate responses.

(11)

Use regression analysis to collect all the effects for each product/process characteristic from step 10.

(12)

Confirm the prediction results by examining the repeatability in the hierarchy of the replicated recipes from step 9.

2.4. The Computational Aids

The basic regression analysis to test replication adequacy, the boxplot response graphs, the correlation analysis for the multivariate responses, the normality tests, as well as the associated QQ plots for the product/process characteristics data, the data shape statistics, and the standard linear and quadratic modeling estimations (curve fitting and ANOVA treatments) were all carried out using the statistical software package IBM SPSS v.29.

To proceed with the specialized computational work, the respective routines were executed on the statistical freeware platform R (v. 4.3.0) [60,61]. The non-linear L₉(3⁴) OA array was constructed using the module ‘param.design()’ from the R-package ‘DoE.base’ (v. 1.2-2). To obtain the Pareto frontier solutions, the R-package ‘rPref’ (v1.4.0) was utilized. The module ‘plot_btg()’ created the ‘Better-than-Graphs’ Hasse diagrams for the preference order of the two-, three-, and four-characteristic response depictions. The module ‘psel()’ was introduced to evaluate the candidate preferences by conducting Skyline queries, along with extending the solution horizon by also obtaining weighted multi-response solutions. The ‘plot_front()’ module provided the two-dimensional skyline graphs. The linear regression via Gibbs sampling was carried out by employing the R-package ‘lrgs’ (0.5.4) [62]. In the module ‘Gibbs.regression()’, the Gibbs sampler was run 1100 times, and the first 100 points (burn-in period) were not included in the simulated conditional posterior distributions of the four product/process characteristic responses.

3. Results

3.1. Data Collection and Screening

The three-level orthogonal recipes of Table 2 were conducted, and the prescribed manufacturing trial runs were duplicated. The responses of the four characteristics (energy consumption, gas consumption, cutting time, and surface roughness) were collected and tabulated in

Table 3. Duplicated L₉(3⁴) OA-designed output for the four product/process characteristics: energy consumption (EC), gas consumption (GC), cutting time (CT), surface roughness (SR).

Run#	EC1	EC2	GC1	GC2	CT1	CT2	SR1	SR2
1	8.44	8.44	7.50	7.50	38	38	5.732	5.8
2	4.89	4.67	19.43	20.35	22	21	5.994	5.992
3	3.78	3.56	33.52	35.62	17	16	6.034	5.932
4	23.75	23.75	11.25	11.25	38	38	5.962	5.956
5	15.00	15.00	23.75	23.75	24	24	6.054	5.952
6	12.50	11.88	14.25	15.00	20	19	6.052	6.066
7	44.33	40.83	15.00	16.28	38	35	6.096	6.02
8	28.00	24.50	11.87	13.57	24	21	4.546	4.542
9	23.33	21.00	21.37	23.75	20	18	5.234	5.612

. In Figure 4A, we have assembled the entire laser-cut specimens. In Figure 4B, we show a close-up side view of two specimens, such that the surface quality is visually discerned by observing the extent and the formed patterns of the striations. In

Table 4. Model summary of the linear regression (standardized) coefficients for the four duplicated characteristic responses (IBM SPSS v.29).

	Standardized Coefficients
Characteristics	βo (p-Value)	β1 (p-Value)	Adj. R2
EC	0.557 (0.361)	0.997 (<0.001)	0.993
GC	0.004 (0.996)	0.996 (<0.001)	0.990
CT	−1.646 (0.295)	0.991 (<0.001)	0.979
SR	0.734 (0.219)	0.962 (<0.001)	0.914

, the model summary of the standardized linear-regression coefficients (IBM SPSS v.29) is listed to assess the practical adequacy of the repeatability of the response data for the four process/product characteristics owing to the trial duplication. All four characteristic intercepts (β_o estimations) are not statistically significant (α = 0.05), thus there is no ‘intercept drift’ between the duplicated runs. Moreover, their slopes (β₁ estimations) are very close to unity for all four process/product characteristics, an indication that there is no ‘sensitivity drift’ between the duplicated runs. Those results are encapsulated by their achieved high fitting performance, as demonstrated by the respective high scores of their adjusted coefficients of determination, which were estimated to be between 0.914 and 0.993. At this point, we may infer that a mere duplication was sufficient to satisfactorily describe the precision of the measurements for all nine OA trial formulations.

However, Figure 5 (IBM SPSS v.29) provides additional insight about the inner trends of the four duplicated response datasets. The first noticeable observation is that the groupings of the two series of trial runs are rather asymmetric, whereas the direction of the asymmetry is always consistent to both trial sets in all four process/product characteristics. This asymmetry appears to be prominent for the gas consumption, the cutting time, and surface roughness characteristics. The 95% confidence intervals for all four plots are depicted to be substantially narrow. This occurrence allows one or two data points to lie outside this range, while several other data points are frequent near the interval borders. Consequently, the inspection checks on this basic data-screening phase suggest the worthwhile engagement of a proper robust profiler to better handle the deviating propensities of the data points by a multi-response multi-parameter optimizer.

The second round of the basic data-screening procedure involves boxplot factorial contrasting, for each type of the four process/product characteristic response datasets, compartmentalized according to their three prescribed factorial settings. It is a quick and practical approach to discern possibly essential predictors from other dormant influences. It is a supplementary assistance to the decision-making process to initiate the more formal multi-parameter statistical treatment of the effects, which follows. Figure 6 (IBM SPSS v.29) gathers all partial plots that quantify robust central tendencies (medians) and variabilities (interquartile ranges) for each process characteristic, and for each sorted controlling-factor setting.

The first remark is that key influences are identified to the individual relationships between the following: (1) the laser power and the energy consumption (Figure 6A), (2) the gas pressure and the gas consumption (Figure 6G), and (3) the cutting speed and the cutting time (Figure 6J). It is noted that from the three mentioned relationships, only the third (CT vs. CS) may be traced to a curvilinear behavior. The other two (EC vs. LP and GC vs. GP) might be simplified by assuming a linear relationship. Another obvious finding is that most data points show strong asymmetry when they are grouped according to their factorial settings. This behavior may be all-embracing for some characteristic responses (EC vs. LP, EC vs. GP, GC vs. F, CT vs. LP, CT vs. F, SR vs. CS, and SR vs. GP), or partially manifested in one-setting (SR vs. LP and SR vs. F) or even two-setting situations (EC vs. F, GC vs. LP, GC vs. GP and, CT vs. GP). The befalling skewness categorization may neither be unidirectional among different settings nor for any given process characteristic. For example, the variety of the median locations is distinct in the subplots C, D, and H of Figure 6. Similarly, another finding is that the magnitude and the spread of the interquartile range (IQR) estimations vary among different factorial settings. There is a broad range of IQR estimate fluctuations that spans from a minimal spread in the CT vs. CS graph, in the first factorial setting (Figure 6J), to include much wider spreads, as exhibited in the graphs of CT vs. LP (Figure 6I), CT vs. GP (Figure 6K), CT vs. F (Figure 6F), SR vs. LP (Figure 6M), SR vs. CS (Figure 6N), SR vs. GP (Figure 6O), and SR vs. F (Figure 6P). Additionally, there is a great variation in the ‘within-a-characteristic’ IQR estimates, which is portrayed with a perceptible anisotropy in such response graphs as EC vs. LP, EC vs. F, GC vs. LP, GC vs. CS, GC vs. GP, GC vs. F, SR vs. LP, SR vs. CS, SR vs. GP, and SR vs. F. However, in the cases of the response graphs of CT vs. LP, CT vs. GP, and CT vs. F, there is a greater similarity in the IQR variability estimations across the profiled factorial settings for each individual process characteristic.

In multi-response screening/optimization studies, correlations between the examined process/product characteristics may be crucial in refining the accuracy of predictions. Customarily, in a correlated group of characteristics, only one characteristic is permitted to enter the data processing phase of the optimizer. Nevertheless, if correlated characteristics are deemed important, because of providing significant variability in the response dataset, then, it may be useful to examine how they may alter the overall factorial-profiling performance; it is especially interesting to understand how they influence the modification of the optimal control settings. In such a situation, they might participate in the multivariate solver arrangement to furnish an augmented “weighted” contribution that serves to assign greater representation to the correlated characteristics. The optimal settings should be determined by the majority of influences that point toward the same direction.

Ordinary correlations have been computed in pairs for the four process/product characteristics. The pertinent correlation estimations are given in 95% confidence intervals in

Table 5. Correlation coefficients between the four process/product characteristics and their 95% confidence intervals (IBM SPSS v.29).

Variable 1	Variable 2	Correlation	Count	Lower C.I.	Upper C.I.
CT	EC	0.456	18	−0.014	0.761
	GC	−0.701	18	−0.880	−0.348
	CT	1.000	18	--	--
	SR	0.184	18	−0.310	0.599
EC	EC	1.000	18	--	--
	GC	−0.349	18	−0.702	0.141
	CT	0.456	18	−0.014	0.761
	SR	−0.238	18	−0.635	0.257
GC	EC	−0.349	18	−0.702	0.141
	GC	1.000	18	--	--
	CT	−0.701	18	−0.880	-0.348
	SR	0.226	18	−0.269	0.627
SR	EC	−0.238	18	−0.635	0.257
	GC	0.226	18	−0.269	0.627
	CT	0.184	18	−0.310	0.599
	SR	1.000	18	--	--

(IBM SPSS v.29). The only significant correlation is identified to the paired association between the gas consumption (GC) and the cutting time (CT); the 95%-upper/lower confidence interval bounds are located between −0.880 and −0.348. However, this range indicates the possibility of this association to be strong as well as weak. Moreover, it is recognized as a negative correlation. This means that as one variable increases, the other one has to decrease.

Incidentally, the generic screening/optimization goal is that both characteristics, the gas consumption and the cutting time have to be synchronously minimized. To resolve this issue, a simple remedy may be to retain both characteristics. That way, a compromised solution may be sought by the multivariate optimizer, which attempts to balance the countering effect of the two characteristic responses.

3.2. Pareto Frontier Analysis

The multi-response multifactorial profiling problem has been treated using a Pareto frontier multi-objective optimizer. Results have been obtained to satisfy two perspectives: (1) to provide a full multi-objective optimal solution, and (2) to partition the solution by including only a partial group of characteristics.

The former objective is achieved by tracking the respective Hasse diagrams, and the latter by plotting 2D Pareto Skylines. The direction of process/product improvement was guided by the minimization goals of the energy and gas consumptions, as well as the minimization of the cutting time, insofar as the surface roughness performance was maintained within a customary laser cutting range of 0.8–6.3 μm. In Figure 7, the Hasse diagrams for the preference ordering (‘Better-than-Graph’ [BTG] version) have been prepared so as to generate two-, three-, and four (full)-objective screening/optimization solution alternatives. All trial runs were input in the module data frame as single records. This tactic was preferred to having the duplicated response dataset split and reduced to central tendencies (means) and variability (signal-to-noise ratio) components, a Taguchi-method basic response-dataset treatment. This action would otherwise double the total number of objectives, and it would inflate the perplexity of the task for the implemented multi-objective solver. All 18 observations were fed into the BTG solver by just four separate response vectors. Thus, it was offered a chance to fine-tune its searching process. Additionally, it introduced a self-confirming step to pinpoint the optimal factorial-setting pattern. This was attained by exposing the paced sequence of the better-performing factorial combinations and their transitive associations between trial runs. The master (4-response) solution may be traced to the Hasse diagram in Figure 7K. It indicates that the optimal input parameters are linked to the first L₉(3⁴) OA trial run (trial #1). Automatically, the solution is confirmed by its duplicated trial outcome (trial #10). Therefore, the emerging final (profiled) adjustments were as follows: (1) the laser power to be set at 4 kW, (2) the cutting speed to be set at 0.8 mm/min, (3) the gas pressure to be set at 0.5 bar, and (4) the focus length to be set at −4 mm. Correspondingly, the achieved performances were as follows: (1) energy consumption of 8.44 kWh, (2) gas consumption of 7.50 Nm³, (3) cutting time of 38 s, and (4) surface roughness estimate at 5.732 μm. It must be remembered that this is a solution in which the correlated gas consumption characteristic was not removed from the screening/optimization scheme, i.e., all four characteristics were present in the BTG layout.

Fortuitously, the final optimal solution remained steadfast regardless of eliminating either the gas consumption or the cutting time. Thus, from Figure 7J, we obtain the same final solution, which is again confirmed by the next in hierarchy better-performing trial-run outcome; it is again tracked to the #10 trial-run. Based on the information which is gleaned from Figure 7G–I, the participation of the gas consumption in the screening scheme merely strengthens the proposed solution against the rest of the ‘three-way’ outcome alternatives. Similarly, with respect to the corresponding ‘two-way’ optimized solutions, the emerging solution is identical to the previous scenarios, according to Figure 7A (EC vs. GC), Figure 7D (GC vs. CT), and Figure 7F (CT vs. SR).

A compromised multi-objective solution for the entire group of examined characteristics was also tested for stability from a different computational perspective. The conditions of achieving parallel minimization goals under the influence of two negatively correlated characteristics were not altered. An expression function created the preferences utility (‘psel()’ function in R-package rPref() v.1.4.0). A weighted scoring was defined for all four individual characteristics which entered the overall minimization goal. The resulting solution pointed to the trial outcome #11, which modifies, in brief form, the minimized characteristic performances as follows: (EC, GC, CT, SR) = (4.67 kWh, 20.35 Nm³, 21 s, 5.992 μm). This new prediction greatly improves the productivity performance, by 81%, at the expense of higher gas consumption. It was also noted that the energy consumption was synchronously reduced by 81%. This improved performance was achieved by setting the inputs (LP, CS, GP, F) to (4 kW, 1.9 mm/min, 0.75 bar, +2.25 mm). The new recommendation agrees only on the laser power optimal setting, in comparison to the previously proposed solutions.

To complete the optimal solution assessment, the portrayals of the Pareto Skylines offer extra information about the limitations of the suggested solutions that we discussed above. For practical purposes, only 2D Pareto frontiers were prepared. All six two-way characteristic combinations are drawn in Figure 8. Following up the Pareto-front prioritized solution, owing to the minimization criteria for the energy consumption and the cutting time, the tradeoff (efficient) solution location appears to be situated in a region (CT vs. GC plot in Figure 8) where its upper bound is approximately less than 15 Nm³ and 20 s, respectively.

From the Pareto Skyline plot, the concurrent minimization of the cutting time and the energy consumption suggests an efficient prediction in the neighborhood of 15 s and about 3 kWh, respectively. If energy and gas consumptions are minimized, the dominating performance would be in the proximity of 10 kWh and 6 Nm³, respectively. The rest of the Pareto front plots aided in ascertaining the tolerated range for the surface roughness performance, which was found to be located below a cutoff value of 6 μm, in all depicted query combinations.

3.3. Gibbs Sampling and Linear Regression Analysis

The saturated multivariate L₉(3⁴) OA dataset underwent a linear regression treatment to approximate the strength of its modelled parameters. The Gibbs sampling method was employed to retrieve the posterior conditional distributions for the four-response multi-parameter problem. The linear relationship between regressors and responses was assumed based on the preliminary visual evidence which was obtained from Figure 6. Therefore, the cause-and-effect relationship could be quantified by a simple linear regression modelling effort, since the regression coefficients could be estimated from the posterior conditional distributions. In

Table 6. Regression coefficients for the concurrent multi-response Gibbs sampler (four process/product characteristics and four controlling factors).

	Process and Product Characteristics
Controlling Factors	EC		GC		CT		SR
	Mean	se	Mean	se	Mean	se	Mean	se
Constant Term	−2.27	0.15	−5.29	0.097	46.89	0.158	6.09	0.02
LP	12.32	0.037	−1.81	0.024	0.29	0.038	−0.29	0.005
CS	−6.06	0.037	6.17	0.025	−9.49	0.038	−0.04	0.005
GP	2.36	0.036	6.57	0.025	−0.52	0.039	0.27	0.005
F	1.35	0.036	0.75	0.025	−0.64	0.037	−0.11	0.005

, the mean coefficients of the four controlling factors are tabulated along with their respective model intercepts for each of the four process/product characteristics. The complete multivariate landscape of the sampled posteriors, expressed as a collection of conditional distributions, was established and used to estimate the regression coefficients (Figure 9).

From Table 6, it may be inferred that the laser power is the predominant effect; its increase induces greater energy consumption. On the other hand, the cutting speed also contributes as the second effect in the strength hierarchy; higher cutting speeds reduce overall energy consumption. Increasing the cutting speed and the gas pressure promotes gas consumption surges at approximately the same rate. It is clear that reducing cutting time is only attained by increasing cutting speed. The quality performance of the surface roughness does not appear to be modulated in the tested input range of any of the examined controlling factors in this study.

4. Discussion

The results that were obtained in the previous section may need additional probing that focuses on the validity of the limitations of the adopting simplifications. In spite of prepping up the multi-objective profiler to handle a non-linear OA mini-dataset, simple linear regression methods were eventually used in conjunction with the Gibbs sampler. Usually, the traditional data analysis of orthogonal datasets commences with the implementation of rudimentary methods that include a data comparative summary in terms of response graphs and response tables. Response graphs and tables are customarily not accompanied with statistical bounds on their estimated quantities; there are only response location estimates on the pre-assigned factorial settings. Consequently, a multifactorial treatment like ANOVA or GLM is necessitated to supply the formal statistical inference of the investigated effects. All of the mentioned approaches rely on the normality of the collected dataset to be evident so as to allow their predictions to become applicable.

The data inspection commences by testing the visual location and dispersion traits of the overall trends of the four process/products characteristics using the QQ-plot and boxplot graphics (Figure 10). From Figure 10A,B (IBM SPSS v.29), it appears that the expected normal behavior for the observed data of the energy and gas consumption are matching fairly well; the boxplot for the gas consumption data may display some minor deviation from symmetry. On the other hand, the cutting time and the surface roughness datasets (Figure 10C,D) exhibit considerable departure from normality in both types of visuals. The boxplot for the cutting time depicts the location estimate strongly trending toward the lower readings; the QQ-plot seems to confirm that behavior. The QQ-plot for the surface roughness response reveals an abrupt shift that is identified to the lower (more sporadic) values of measurements. The boxplot for the surface roughness makes vivid the data asymmetry in connection to higher values that otherwise are marked by small variability.

A check on the normality status of the four process/product characteristics is shown in

Table 7. Normality tests for the four process/product characteristics: energy consumption (EC), gas consumption (GC), cutting time (CT), and surface roughness (IBM SPSS v.29).

Characteristics	Kolmogorov–Smirnov a			Shapiro–Wilk
	Statistic	df	Sig.	Statistic	df	Sig.
EC	0.142	18	0.200 *	0.909	18	0.083
GC	0.150	18	0.200 *	0.925	18	0.157
CT	0.267	18	0.001	0.803	18	0.002
SR	0.309	18	<0.001	0.662	18	<0.001

* This is a lower bound of the true significance. ^a Lilliefors Significance Correction.

. Both standard tests, Kolmogorov–Smirnov and Shapiro–Wilk (IBM SPSS v.29), have been utilized to quantify the level of significance of the results (α = 0.05).

There is agreement on the data normality of the energy and the gas consumption data, whereas this may not be true for the remaining two characteristics. This might not be a serious issue on balancing the contributions among the four characteristics, because the surface roughness estimations were restricted within the acceptable production range. Nevertheless, this deviation tangles the profiler process, possibly ushering in more uncertainty elements in the predictions. It is dubious how the screening results might be interpreted, even when the pivotal estimate of the constant coefficient (intercept) in the model remains computationally questionable. Resorting to a methodology that combined the Gibbs sampling technique with the Pareto frontier multi-objective optimizer ensured that the robust and conditionally estimated regression statistics could overcome any data oddities.

To better understand the peculiarity of the dataset behavior pertaining to the characteristics of the cutting time and the surface roughness, the shape estimators of skewness and excess kurtosis have been computed in

Table 8. Measure of shape statistics (skewness and excess kurtosis) for the four process/product characteristics: energy consumption (EC), gas consumption (GC), cutting time (CT), and surface roughness (IBM SPSS v.29).

	N	Skewness		Kurtosis
	Statistic	Statistic	Std. Error	Statistic	Std. Error
EC	18	0.820	0.536	0.152	1.038
GC	18	0.830	0.536	0.322	1.038
CT	18	0.570	0.536	−1.522	1.038
SR	18	−1.986	0.536	2.981	1.038

(IBM SPSS v.29). The skewness estimates for the energy and gas consumptions, as well as the cutting time, basically imply symmetric distribution forms, since they are less than a value of one. However, given the fact that the standard error has been computed to a value of 0.536, the right-skewness of the dataset cannot be rejected. On the other hand, the surface roughness estimate indicates that it is more likely to follow a left-skewed distribution. Excess kurtosis for energy and gas consumptions is reasonably small (close to normal) for both characteristics. Still, their larger standard error values might not exclude other possible manifestations (platykurtic/leptokurtic trends). Conversely, the excess kurtosis estimates for the cutting time and the surface roughness are more dichotomized, the former favoring a more platykurtic shape and the latter a more leptokurtic shape. Again, the other possible classifications cannot be ruled out owing to their large standard error estimation values.

It is informative to discuss the validity of the assumption to simplify the impending non-linear optimizer to an approximate linear model. The ensuing statistical analysis is organized by separating the multifactorial treatment for each individual process/product characteristic. A comparison of the statistical estimations of the regression coefficients, their p-values, the accompanying ANOVA statistical significances, and their adjusted coefficients of determinations has been separately arranged in

Table 9. Linear and quadratic modeling of the four process/product characteristics for each of the four controlling factors (IBM SPSS v.29): Statistical significances (p-values) for ANOVA and regression coefficients and their corresponding adjusted coefficient of determination (R²).

		Linear Model			Quadratic Model
Characteristic	Factor	ANOVA sig.	Coefficient sig.	Adjusted R2	ANOVA sig.	Coefficient sig.	Adjusted R2
EC	LP	<0.001	<0.001	0.72	<0.001	0.541	0.703
						0.765
	CS	0.078	0.078	0.13	0.187	0.405	0.094
						0.557
	GP	0.509	0.509	−0.033	0.798	0.936	−0.1
						0.862
	F	0.713	0.713	−0.053	0.819	0.579	−0.104
						0.611
GC	LP	0.439	0.439	−0.022	0.643	0.52	−0.069
						0.588
	CS	0.003	0.003	0.398	0.013	0.421	0.363
						0.734
	GP	0.002	0.002	0.442	0.008	0.658	0.405
						0.94
	F	0.742	0.742	−0.055	0.847	0.675	−0.108
						0.638
CT	LP	0.897	0.897	−0.061	0.939	0.734	−0.124
						0.745
	CS	<0.001	<0.001	0.878	<0.001	<0.001	0.971
						<0.001
	GP	0.847	0.847	−0.06	0.982	0.979	−0.131
						1
	F	0.796	0.796	−0.058	0.962	0.886	−0.128
						0.914
SR	LP	0.037	0.037	0.196	0.026	0.155	0.304
						0.082
	CS	0.718	0.718	−0.054	0.329	0.144	0.023
						0.154
	GP	0.043	0.043	0.184	0.134	0.591	0.133
						0.808
	F	0.419	0.419	−0.019	0.156	0.109	0.115
						0.084

for the two pertinent versions: the linear model and the quadratic model (IBM SPSS v.29). The first finding is that from the 16 completed regression fittings that were attempted in this examination, all relationships may be fitted to a simple linear expression. It is only the relationship of the cutting time and the cutting speed that may be represented to also a quadratic model. In this case, both linear and quadratic models generate regression coefficients where their p-values are smaller than 0.001. The improvement in explaining the uncertainty appears to favor the quadratic model over the linear version. The estimations of the adjusted coefficient of determination were found to be 0.88 and 0.97 for the linear and the quadratic models, respectively. In practice, the quadratic model could be reducible to the linear model, since the quadratic model predicts a monotonous curvilinear trace (Figure 11D); the optimal setting is expected to be located on the lowest response estimation (largest cutting speed setting). The linear statistical modeling of the energy consumption against the laser power is strongly significant for the linear fitting, but not for the quadratic extension (Table 9). From Figure 11A, it becomes clear that it is superfluous to retain a second order term in the model, which in this case led to a nonsignificant quadratic prediction. Furthermore, from Table 9, it is observed that the linear modeling of the relationships between the gas consumption and the cutting speed, and respectively, between the gas consumption and the gas pressure, return in both cases statistically significant coefficient estimations. In support to this remark, a quadratic model could not justify a distinct relationship, as witnessed from Figure 11B,C. The same findings could be extended to the predicted relationship (p < 0.05) between the surface roughness and the gas pressure (Figure 11F). The linear trend between the surface roughness and the laser power in Figure 11E is still preferable, because its coefficient estimation is statistically significant, while the quadratic approximation returns no statistically significant coefficients at all (p > 0.05). From a practical perspective, the laser power may be adjusted at the lowest operating point (economical solution), which also permits the surface roughness observations to be controlled and be bound to the acceptable production range, which is congruent to the laser cutting capabilities for this type of working material.

5. Conclusions

The study examined the consistency of the surface roughness performance of a carbon EN10025 St37-2 10mm-thick workpiece that was processed by a Trumpf 5040 CO₂ 2D-TrueLaser cutting unit. The surface roughness performance was fairly repeatable and well within the typical range of 0.8–6.3 μm which is expected in the metals processing industry; the conducted trial measurements varied from 4.542 to 6.096 μm. However, the caveat of this application stemmed from the parallel benefits of inculcating lean-and-green principles to the metal processing task. The novelty of this work lies in demonstrating how to address two distinctly different opportunities for reaping lean-and-green savings, as part of the product/process improvement. One opportunity relates to cost reduction, emanating from implementing an efficient experimental tactic, and the other on the lean-and-green datacentric multi-objective process/product development. The first exploited opportunity regards the dramatically reduced demand for resources (materials and equipment availability), which ensured the feasibility of this project. Implementing a non-linear Taguchi-type L₉(3⁴) OA minimized the total logistics requirements down to only 1/9 of the scheduled trials for an otherwise full-factorial experimental scheme. The motivation to resort to such a lean-and-green tactic is well-aligned with the UN Sustainability Development Goals #9 and 12.

The emphasis for this process/product screening/optimization exercise was placed on synchronously minimizing energy and gas consumptions, while boosting the productivity performance (lower cutting time), besides achieving and maintaining high repeatability in the primary surface quality-integrity characteristic. This required data collection for the energy and gas consumptions on the laser cutting machine center, as well as the timing for each complete cutting unit cycle. The introduction of robust tools to handle the resulting multi-response problem was pivotal as Taguchi-type orthogonal samplers generated the small-and-dense datasets. The resulting data properties may or may not be comfortably treatable with ordinary techniques. This occurs perhaps for two reasons. First, any data concerns about their distribution details usually remain unresolved. This is due to the high uncertainty which is congenital to situations when one is confined to work with small samples. Second, the inherent high complexity issue which underlies the synchronous manipulation of multiple effects to control the performance of diverse characteristics. Consequently, a double arrangement of multi-objective solvers, the Pareto frontier approach and the Gibbs-sampling regression analysis, was attempted in the ensuing multivariate multifactorial screening/optimization scheme in order to ensure that the strength of the four investigated controlling factors (laser power, gas pressure, cutting speed, and focus length) were statistically determined. The Pareto frontier optimal solution was (EC, GC, CT, SR) = (4.67 kWh, 20.35 Nm³, 21 s, 5.992 μm) for the four-characteristic minimization objective. This corresponds to factorial settings adjusted at the inputs (LP, CS, GP, F) to (4 kW, 1.9 mm/min, 0.75 bar, +2.25 mm). This solution was strengthened by the multi-response posterior distributions of the Gibbs regression analysis that associated (1) the positive effect of the laser power and the negative effect of the cutting speed, on energy consumption, (2) the positive effect of the cutting speed and the gas pressure on gas consumption and, (3) the negative effect of cutting speed on cutting time. It is an interesting fact and a novel occurrence to this multivariate profiling that the multi-response solver scheme retained and evaluated two negatively correlated characteristics, the gas consumption and the cutting time. This condition had to elicit a compromised multivariate solution.

The minimization of the energy and gas consumptions, along with the maximization of productivity, called for an innovative industrial solution (UN SDG #9). The optimization of the laser cutting performance properties “ensures sustainable consumption and production patterns” by ensuring good use of resources which also improves energy efficiency (UN SDG #12). Further work can be easily extended by exploring different materials, even complex composites, along with a variety of test unit dimensions in an assortment of various geometrical shapes.