Evolutionary Optimizing Process Parameters in the Induction Hardening of Rack Bar by Response Surface Methodology and Desirability Function Approach under Industrial Conditions

Conditions of industrial production introduce additional complexities while attempting to solve optimization problems of material technology processes. The complexity of the physics of such processes and the uncertainties arising from the natural variability of material parameters and the occurrence of disturbances make modeling based on first principles and modern computational methods difficult and even impossible. In particular, this applies to designing material processes considering their quality criteria. This paper shows the optimization of the rack bar induction hardening operation using the response surface methodology approach and the desirability function. The industrial conditions impose additional constraints on time, cost and implementation of experimental plans, so constructing empirical models is more complicated than in laboratory conditions. The empirical models of nine system responses were identified and used to construct a desirability function using expert knowledge to describe the quality requirements of the hardening operation. An analysis of the hypersurface of the desirability function is presented, and the impossibility of using classical gradient algorithms during optimization is empirically established. An evolutionary strategy in the form of a floating-point encoded genetic algorithm was used, which exhibits a non-zero probability of obtaining a global extremum and is a gradient-free method. Confirmation experiments show the improvement of the process quality using introduced measures.


Introduction
Engineering requirements of gear parts that have specific properties lead to research on developing the induction hardening processes.Especially in the automotive industry, reaching a compromise between cost, quality and reliability by using semi-automatic heat treatment processes to ensure the mechanical and geometrical features is receiving more and more attraction.However, the induction hardening process is challenging due to thermal strains affecting the functionality and quality of crucial elements in the produced structures.Other limitation factors are the formation of non-hardened zones in parts and the change in the designed hardness of the treated elements.Moreover, the induction hardening process is quite complex; it could be represented as a coupled field problem with mechanical, electromagnetic and thermal fields, which are pretty hard to model in the industrial environment, especially in the serial mode of production.
Optimization of the induction hardening process is well-known in the literature.Two main approaches can be distinguished here: one using complex computational mechanics models and the other using empirical modeling.The following works represent the computational mechanics approach.Favennec et al. [1] show how to solve the problem of temperature distribution optimization using optimal control techniques, and similarly, Jakubovicowa et al. [2] optimized the process for the uniform surface temperature distribution criterion.Nemkov et al. [3] simulated the stress and distortion evolution during the induction hardening process of tubes using the finite element method; the hardness profile sensitivity in the induction hardening process with the finite element simulations was explored in [4].The coupled electromagnetic, thermal and mechanical computational models of the heat treatment process using electromagnetic fields were shown in [5].Fisk et al. [6] have presented the complex models of induction hardening in low alloy steels.Reference [7] shows the optimization of the edge effect of 4340 steel specimens heated by induction process with flux concentrators using finite element axis-symmetric simulation.The computational approach brings several advantages, such as universality and flexibility, as shown in [8,9], and also lower cost and faster results than an experimental and analytical approach.However, the computational approach requires tuned, complex constitutive models, which are hard to identify, especially in industrial conditions.That is why the empirical modeling approach was also successfully applied to solving various optimization problems [10], also arising during the induction hardening process.
The empirical modeling approach allows for establishing relations between process outputs and inputs using experimentation and measurements with statistical techniques for managing data [11].Kohli and Singh [12] have used response surface methodology (RSM) to find the optimal values of process parameters for induction hardening of AISI 1040 steel.Various process parameters, such as feed rate, current, dwell time, and the gap between the workpiece and induction coil, are experimentally explored.In [13], the RSM and Taguchi method optimized the induction hardening process for maximum depth and minimum edge effect.The multiobjective optimization problem with the appropriate economic, environmental, and social metrics was analyzed in [14] to assure the sustainability of the induction hardening process using empirical models.The reduction of edge effect using the RSM and artificial neural network modeling of a spur gear treated by induction with flux concentrators was shown in [15].Artificial intelligence modeling of induction contour hardening of 300M steel bar and C45 steel spur-gear was performed in [16].Reference [17] shows that the central composite design, with a second-order response surface design, was employed to systematically estimate the empirical models of temperature and phase transformation geometry during the induction hardening.The effect of scanning speed and air gap on the uniformity of hardened depth and mechanical properties of large-size spur gears was investigated in [18].Multi-response optimization using the desirability function approach of the induction hardening process using quality responses such as the effective case depth and hardness values were analyzed in [19] for different combinations of medium frequency power, feed rate, quench pressure and temperature.Asadzadeh et al. [20] have shown the hybrid model, integrating measurements and physics, of the induction hardening.
A gap in the literature can be identified based on the authors' best knowledge and the review presented.It concerns the solution to the problem of multi-criteria optimization of the induction hardening process in industrial conditions utilizing a hybrid approach using empirical modeling and computational intelligence tools.The novelty of this work lies in the application of multiple qualitative metrics regarding deformation, hardness profile and hardening depth in a complex structural component such as a steering gear rack bar to optimize the parameters of the induction hardening process, using empirical modeling, a desirability function and an evolutionary algorithm.
Some essential points characterize the present work, which are identified and presented in this paper:

•
There is no evidence that evolutionary optimization using a genetic algorithm has been applied to induction hardening processes under complex industrial process requirements, modeled empirically; • Lack of availability of studies showing the effect of the form of the desirability hypersurface on the effectiveness of optimization algorithms for the induction hardening process; • Lack of availability of studies showing the effectiveness of global optimization algorithms, such as genetic algorithms, for obtaining a set of quasi-optimal solutions to the problem for the desirability function formulated for the induction hardening process;

•
There is a lack of availability of studies showing the complex problem of conducting experiments to confirm the optimal parameter settings of the induction hardening process under industrial conditions and analyzing their results contained in small samples, which precludes the use of parametric tests.
In this work, RSM with the central composite design (CCD) of experiments was employed to establish the functional relationship between the three main process parameters that served as design variables of the induction hardening multi-criteria optimization problem.Several responses related to hardness profiles and geometrical measures of quenching quality of the heat treatment operation of the automotive steering gear were formulated.Among them, the most important in the following work is the problem of minimizing residual thermal deformations, which is enforced by an additional straightening operation, increasing the duration of the process and its costs.However, the desirability approach allowed for transforming the multiobjective optimization problem into a single-optimization problem using expert knowledge for setting the weights in the desirability function.The quadratic models for the process responses were quantitatively analyzed, and their significance and accuracy were confirmed statistically.Next, the reduction of the models is performed to consider significant terms of the regression models.This step introduces a specific consequence for the optimization problem formulation in a change of an optimized function form.That change affects the effectiveness of the applied optimization algorithms.The article shows that the global optimization technique as the evolutionary strategy in the form of GA allows to meet the difficulties resulting from the form of the optimized function.Finally, the confirmation experiments should be conducted to verify the optimal solution the GA method determines.
This article consists of six Sections.After the introduction in the present Section, Section 2 briefly describes the considered induction hardening process with its quality indicators and experimentation methods.Section 3 is devoted to a description of the empirical models that were obtained with RSM.Section 4 focuses on formulating the multi-criteria optimization problem and obtaining the solution using the computational intelligence technique, namely the GA algorithm.Section 5 analyzes confirmation experiments.The last Section, Section 6, briefly summarizes the conducted research.

Induction Hardening Process
The process is conducted using a particular equipment.The automatic inductive hardening and tempering machine have a rotation table with three stations, as shown in Figure 1:

•
For loading rack bars before the heat treatment process and unloading rack bars after the process is finished; The machine is available for hardening and tempering rack bars with lengths of 500-900 mm and 22-32 mm diameters.
The hardening machine provides a possibility to control the hardening and tempering process condition by adjusting:

•
The hardening power is in the range of 0-100%;

•
The hardening feed rate ranges from 100-50,000 mm/min;

•
The distance of the hardening coil to rack bar teeth is in the 1.5-5 mm range.The machine is available for hardening and tempering rack bars with lengths of 900 mm and 22-32 mm diameters.
The hardening machine provides a possibility to control the hardening and tem ing process condition by adjusting:

•
The hardening power is in the range of 0-100%;

•
The hardening feed rate ranges from 100-50,000 mm/min; • The distance of the hardening coil to rack bar teeth is in the 1.5-5 mm range.
Moreover, the flow of quenching water could be controlled; however, due to th accuracy of the machine flow meter and lack of control possibility in auto mod quenching flow rate was kept constant through the process with a value equal to 45 L Moreover, the hardening station has sensors in order to provide a continuous mea ment of a rack bar distance from the coil.It also provides a stable distance from the ening coil during the hardening process.Furthermore, during the tempering pro which gave the required surface hardness, the power and the feed were kept con through experiments with settings:

•
The feed rate equals 1600 mm/min.
The measurements of the following quantities were conducted during the qu control of rack bars: 1. Hardness depth; 2. Surface hardness; 3. It is mandatory to use a straightening process after hardening due to deform which is an effect of the hardening process; the time of the straightening oper increases when rack bars deformation is over 1100 µm.Moreover, the flow of quenching water could be controlled; however, due to the low accuracy of the machine flow meter and lack of control possibility in auto mode, the quenching flow rate was kept constant through the process with a value equal to 45 L/min.Moreover, the hardening station has sensors in order to provide a continuous measurement of a rack bar distance from the coil.It also provides a stable distance from the hardening coil during the hardening process.Furthermore, during the tempering process, which gave the required surface hardness, the power and the feed were kept constant through experiments with settings:

•
The feed rate equals 1600 mm/min.
The measurements of the following quantities were conducted during the quality control of rack bars:
It is mandatory to use a straightening process after hardening due to deformation, which is an effect of the hardening process; the time of the straightening operation increases when rack bars deformation is over 1100 µm.
The main aim of the present study was to establish hardening process parameters in order to achieve hardening conditions that match quality requirements, i.e., minimize rack bar deformations with the proper hardening depth and hardness profiles and avoid increasing a total machine cycle time, mainly when it includes additional operations, e.g., straightening.
Additional ancillary factors, e.g., coolant state, coolant flow rate, environment temperature, humidity, hardening coil condition, deviation in allow composition of steel, etc., can influence the output variables; nevertheless, the performed experiments are restricted to the three hardening parameters, which are automatically controlled by the machine.Figure 2 shows the process setup.
Materials 2023, 16, x FOR PEER REVIEW 5 of 24 increasing a total machine cycle time, mainly when it includes additional operations, e.g., straightening.Additional ancillary factors, e.g., coolant state, coolant flow rate, environment temperature, humidity, hardening coil condition, deviation in allow composition of steel, etc., can influence the output variables; nevertheless, the performed experiments are restricted to the three hardening parameters, which are automatically controlled by the machine.Figure 2 shows the process setup.The hardening operation in automatic mode is as follows.After the hardening element is delivered to the loading station, it is mounted in holders, and the inductor is fixed in the initial position at the bottom, at a predetermined distance from the hardening piece, according to the plan of experiments.At a fixed power and feed rate, the quenching process begins, during which the coil moves upwards towards the upper holder with the coolant flow, as shown in Figure 2.After passing the set distance, the holder is released, the coil is discharged to the neutral position, and the hardened element is transferred to the tempering operation.

Process Quality Indicators and the Measurement Setup
Specific measurements are required to identify the quality of the hardening process qualitatively.Below, the list of performed measurements with the quality requirements is presented:

•
The hardness depth on the teeth side with minimum requirements equal to 3.9 mm; • The hardness depth on the back side (an opposite side to teeth) with requirements given by the range above 1 mm;

•
The surface hardness on teeth with requirements given by the range 55-60 HRC;

•
The surface hardness on the back side is within the 52-55 HRC range requirements; • After finishing the hardening and tempering operations, all rack bars are straightened to obtain a maximal deflection of value not greater than 1000 µm before going to the next operation; in the presented study, the existence of thermal strains after induction hardening and tempering operations is a primary driving force for performing process optimization to reduce costs and time.
The general flow of the experimentation procedure is presented in Figure 3.
The hardening operation in automatic mode is as follows.After the hardening e ment is delivered to the loading station, it is mounted in holders, and the inductor is fix in the initial position at the bottom, at a predetermined distance from the hardening pie according to the plan of experiments.At a fixed power and feed rate, the quenching p cess begins, during which the coil moves upwards towards the upper holder with t coolant flow, as shown in Figure 2.After passing the set distance, the holder is releas the coil is discharged to the neutral position, and the hardened element is transferred the tempering operation.

Process Quality Indicators and the Measurement Setup
Specific measurements are required to identify the quality of the hardening proc qualitatively.Below, the list of performed measurements with the quality requirements presented:

•
The hardness depth on the teeth side with minimum requirements equal to 3.9 mm The hardness depth on the back side (an opposite side to teeth) with requireme given by the range above 1 mm; The surface hardness on teeth with requirements given by the range 55-60 HRC; The surface hardness on the back side is within the 52-55 HRC range requiremen • After finishing the hardening and tempering operations, all rack bars are straig ened to obtain a maximal deflection of value not greater than 1000 µm before goi to the next operation; in the presented study, the existence of thermal strains af induction hardening and tempering operations is a primary driving force for p forming process optimization to reduce costs and time.
The general flow of the experimentation procedure is presented in Figure 3.The zones associated with the hardness measurements are shown in Figure 4. Figu 4 shows six hardness measurement zones: three zones for the tooth side and three for t back side of the rack bar.The numbers indicate which teeth the hardness measureme apply to.Places of the hardness depth measurements are also presented on the cross-secti of the rack bar in the middle of its length, as shown in Figure 5.The zones associated with the hardness measurements are shown in Figure 4. Figure 4 shows six hardness measurement zones: three zones for the tooth side and three for the back side of the rack bar.The numbers indicate which teeth the hardness measurements apply to.The hardening operation in automatic mode is as follows.After the hardening element is delivered to the loading station, it is mounted in holders, and the inductor is fixed in the initial position at the bottom, at a predetermined distance from the hardening piece, according to the plan of experiments.At a fixed power and feed rate, the quenching process begins, during which the coil moves upwards towards the upper holder with the coolant flow, as shown in Figure 2.After passing the set distance, the holder is released, the coil is discharged to the neutral position, and the hardened element is transferred to the tempering operation.

Process Quality Indicators and the Measurement Setup
Specific measurements are required to identify the quality of the hardening process qualitatively.Below, the list of performed measurements with the quality requirements is presented:

•
The hardness depth on the teeth side with minimum requirements equal to 3.9 mm; • The hardness depth on the back side (an opposite side to teeth) with requirements given by the range above 1 mm; The surface hardness on teeth with requirements given by the range 55-60 HRC; The surface hardness on the back side is within the 52-55 HRC range requirements; • After finishing the hardening and tempering operations, all rack bars are straightened to obtain a maximal deflection of value not greater than 1000 µm before going to the next operation; in the presented study, the existence of thermal strains after induction hardening and tempering operations is a primary driving force for performing process optimization to reduce costs and time.
The general flow of the experimentation procedure is presented in Figure 3.The zones associated with the hardness measurements are shown in Figure 4. Figure 4 shows six hardness measurement zones: three zones for the tooth side and three for the back side of the rack bar.The numbers indicate which teeth the hardness measurements apply to.Places of the hardness depth measurements are also presented on the cross-section of the rack bar in the middle of its length, as shown in Figure 5. Places of the hardness depth measurements are also presented on the cross-section of the rack bar in the middle of its length, as shown in Figure 5.
To measure the hardening depth, the micro Vickers method was used, and the FM-810 Micro Vickers Hardness Tester performed the inspection with an accuracy of ±1 HV.This measurement involves determining the microhardness profile along the rack bar's thickness for a tooth cross-section.The measurements assume that the hardening depth defines a point with a hardness of 400 HV.An example of the hardness profile is shown in Figure 6.To measure the hardening depth, the micro Vickers method was used, a 810 Micro Vickers Hardness Tester performed the inspection with an accurac This measurement involves determining the microhardness profile along th thickness for a tooth cross-section.The measurements assume that the hard defines a point with a hardness of 400 HV.An example of the hardness profi in Figure 6.To measure the hardening depth, the micro Vickers method was used, and the 810 Micro Vickers Hardness Tester performed the inspection with an accuracy of ±1 This measurement involves determining the microhardness profile along the rack thickness for a tooth cross-section.The measurements assume that the hardening d defines a point with a hardness of 400 HV.An example of the hardness profile is sh in Figure 6.The Zwick/Roell ZHR 8150LK inspects surface hardness in the HRC scale wi accuracy of ±1 HRC.The mentioned methods are indirect; moreover, destroying the ponent (a rack bar) to perform measurements is necessary.
Runout inspection and straightening are performed by the automatic process rack bar straightening; the "Rack Bar Bend M/C" machine performs measuremen eight of the high accuracy digital displacement transducer gauge probe with an accu The Zwick/Roell ZHR 8150LK inspects surface hardness in the HRC scale with an accuracy of ±1 HRC.The mentioned methods are indirect; moreover, destroying the component (a rack bar) to perform measurements is necessary.
Runout inspection and straightening are performed by the automatic process of a rack bar straightening; the "Rack Bar Bend M/C" machine performs measurements by eight of the high accuracy digital displacement transducer gauge probe with an accuracy of ±1.2 µm.Figure 7 shows how deformation is measured and how the machine performs straightening during the automatic process.

Experimental Design
The central composite design (CCD) [11] as a plan of experiments is used.Three input factors (design variables) were set within the following ranges according to the knowledge of the process before optimization:
Due to the destructive inspection, it is necessary to scrap all components used in the experiments; therefore, due to cost reduction, there is no replication in the current plan of experiments.The results of the experiments are presented in Table 1.Response y1 denotes the maximal deflection of the rack bar after the induction hardening operation is finished.Responses y2 and y3 denote the hardness depths on the teeth and back-side, along the cross-section of the rack bar, in the middle of its length, respec-

Experimental Design
The central composite design (CCD) [11] as a plan of experiments is used.Three input factors (design variables) were set within the following ranges according to the knowledge of the process before optimization:
Due to the destructive inspection, it is necessary to scrap all components used in the experiments; therefore, due to cost reduction, there is no replication in the current plan of experiments.The results of the experiments are presented in Table 1.Response y 1 denotes the maximal deflection of the rack bar after the induction hardening operation is finished.Responses y 2 and y 3 denote the hardness depths on the teeth and back-side, along the cross-section of the rack bar, in the middle of its length, respectively.Responses y 4 -y 6 indicate the average hardness measured in the I, II and III zones on the teeth side and responses y 7 -y 9 indicate the hardness measured in the I, II and III zones on the back side of the rack bar.
Grubb's test [11] was performed for data given in Table 1 to identify outliers.The significance level was 0.05, and the following outliers are identified in the collected data:

•
In experiment no.1-in responses 2, 5 and 6, the result indicates that the deep hardening was obtained without the required hardness of the teeth in the 2nd and 3rd zones; • In experiment no.13-response 4 indicates that the hardness of the teeth in the 1st zone is insufficient; • In experiment no.16-in responses 7, 8 and 9, the result indicates that the hardening of the back side of the rack bar was not achieved; moreover, in response 3-the hardening depth on the back side-is equal to 0, however, on the assumed significance level this response is not an outlier.
It is well known that outliers affect the response surface models, but in the presented study, no replications of experiments were assumed; hence all results were applied for modeling.

RSM for Empirical Modelling
Empirical models of process responses are assumed in the form of full quadratic models as follows: where y (i) is the i-th response, b 0 is the free term of the model, vector x denotes the design variables, vector b collects the linear terms coefficients, and matrix B includes the quadratic terms of the response surface model; it is assumed that the responses are uncorrelated and each model is under the normal noise with 0 mean and standard deviation σ (i) .Models are linear with respect to the coefficients, which is why the least square method (LSM) [11] was applied to identify coefficients.An analysis of variance (ANOVA) is carried out to estimate the validity of the identified mathematical models and the effect of every model term on the responses.The p-value is used to check the significance, which means the response is greatly determined by the model term whose p-value is sufficiently low (less than or equal to the significance level of 0.05).On the contrary, the higher the F-value, the stronger the significance of the model item.Symbol S indicates the significant term in the considered model.
The ANOVA of the mathematical model for the maximal deflection is shown in Table 2.The F-value (5.26) and p-value (0.008) imply that the model is significant, while the F-value (22.24) and p-value (0.002) imply that the lack of fit is significant.If the lack of fit is significant, the higher-order model terms should be taken into consideration to improve the accuracy of the model; however, in the present study, only the full quadratic models are considered.According to the F-value and p-value, there are four significant model terms, among which the power x 1 and the feed rate x 3 have the most significant influence, followed by the linear-by-linear interaction effect between the power x 1 and the feed rate x 3 and the quadratic effect of the power x 1 2 .The coil distance x 2 and other model terms also containing it seem to have little influence on the maximal deflection of the rack bar.The ANOVA of the mathematical model for the hardening depth of the teeth is shown in Table 3.The F-value and p-value of the lack of fit are 0.41 and 0.82, respectively, while the F-value and p-value of the model are 14.41 and less than 0.00013, respectively, which indicates that the model is significant enough and no higher-order model terms need to be considered.The power and the feed rate amount have the most significant influence on the hardening depth of the teeth, followed by the linear-by-linear interaction effect between the power and the feed rate, the coil distance x 2 and the quadratic effect of the power x 1 2 .The linear-by-linear interaction effect between the power and the coil distance is also significant.The ANOVA of the mathematical model for the hardening depth on the back side is shown in Table 4.As can be seen from the table, the model is significant with an F-value (11.71) and a small p-value (<0.0004), while the lack of fit is still not significant with a smaller F-value (0.25) and larger p-value (0.92).The power x 1 has the most significant influence on the hardening depth on the back side, with the largest F-value (69.29), followed by the feed rate x 3 and the linear-by-linear interaction effect between the power and the feed rate.Other model terms seem to have little influence on the hardening depth on the back side of the rack bar.The ANOVA of the mathematical model for the hardness of the teeth in the first zone is shown in Table 5.As can be seen from the table, the model is significant with an F-value of 20.72 and an extremely small p-value <0.0001, while the lack of fit is still not significant with an F-value of 4.42 and a p-value of 0.06.The power x 1 has the most significant influence on the hardness of the teeth in the first zone, with the largest F-value (112.54),followed by the quadratic effect of the power x 1 2 and the feed rate x 3 .The linear-by-linear interaction effect between the power and the feed rate is also significant at the 0.05 significance level.The ANOVA of the mathematical model for the hardness of the teeth in the second zone is shown in Table 6.The model is significant, with an F-value of 10.49 and a p-value < 0.0006, and the lack of fit is still insignificant.The most important factors are the power x 1 , the feed rate x 3 and the quadratic effect of the power x 1 2 .Similar significance levels show the coil distance x 2 and the linear-by-linear interaction effect between the power and the coil distance.The ANOVA of the mathematical model for the hardness of the teeth in the third zone is shown in Table 7.The model is significant, with an F-value of 9.44 and a p-value < 0.0009, and the lack of fit is still insignificant.The most important factors are the power x 1 , the feed rate x 3 and the coil distance x 2 .Similar levels of significance, but with smaller values of the F-statistics than earlier considered factors, show the quadratic effects of the power and the feed rate.The ANOVA of the mathematical model for the hardness of the back-side in the first zone is shown in Table 8.The model is significant with an F-value of 11.32 and p-value < 0.0004, and the lack of fit is significant with a p-value of 0.003.Only three factors are significant: the quadratic effect of the power x 1 , the linear-by-linear interaction effect between the power and the feed rate and the pure quadratic effect of the power.The ANOVA of the mathematical model for the hardness of the back-side in the second zone is shown in Table 9.The model is significant with an F-value of 13.27 and a p-value < 0.0002, and the lack of fit is very significant with a p-value < 0.00002.Only four factors are significant: the quadratic effect of the power x 1 2 followed by the power itself, the linear-by-linear interaction effect between the power and the feed rate, and the feed rate itself.Similar conclusions could be formulated for the mathematical model of the hardness of the back-side in the third zone, as shown in Table 10.The model is significant with an F-value of 14.21 and a p-value < 0.0002, and the lack of fit is significant with a p-value < 0.0007.Moreover, only four factors are significant: the quadratic effect of the power x 1 2 followed by the power itself, the linear-by-linear interaction effect between the power and the feed rate, followed by the feed rate itself.The lack of fit in the presented models indicates the need to use the higher-order models; however, in the present work, only quadratic models will be used.The higherorder models can be applied, e.g., in polynomials, neural networks and kriging.
The identified models are nonlinear, indicating room for improvement or even process optimization.Table 11 presents the essential statistical characteristics of the identified full quadratic models (1).The significance level was assumed as 0.05.The analysis shows that the identified models are strong and significant but with moderate predictive ability, as shown by the PRESS statistics.The backward elimination of the insignificant terms is performed to increase the predictive properties of the identified models.The reduced models will be used in the optimization process; models presented in Table 12 were shown with an accuracy of two significant digits.The values of the PRESS statistics are significantly smaller than for the full models, with a tiny drop in R 2 values.It is essential to emphasize that the reduced models change the background for optimization.Supposedly, we will consider the model for the maximal deflection y 1 .In that case, the canonical analysis shows that for the full quadratic model, the stationary point  12, has a stationary point x s = [−1.0320 −0.377]T which is also a saddle point, but from the ridge system because the eigenvalues are now {−63.10,0, 2246.76}.Figure 8 shows contour plots of the response surface slices at x 3 = 0.5 for the full quadratic model (left) and the reduced model (right).

2
Hardening depth of teeth, y2 It is essential to emphasize that the reduced models change the backgro mization.Supposedly, we will consider the model for the maximal deflecti case, the canonical analysis shows that for the full quadratic model, the sta xs = [−0.1960.764 1.053] T is a saddle point because the eigenvalues of matr 36.93,269.55} are mixed in sign.The reduced model, according to Table 12, ha point xs = [−1.0320 −0.377]T which is also a saddle point, but from the ridge sy the eigenvalues are now {−63.10,0, 2246.76}.The change in the type of the response surface affects the optimization process, especially in the single-objective case.For the multiobjective case, the desirability approach allows taking into account various forms of response surfaces.However, the properties of the ridge systems will also be present in the objective function, as shown in the next Section.

Multiobjective Optimization Problem Formulation and Its Solution
The design variable vector x is as follows x = [x 1 x 2 x 3 ] T , where x 1 is the power, x 2 is the distance between the coil and the part, and x 3 is the feed rate.Then, the multiobjective optimization problem can be formulated using the desirability function D as given below: ( The sphere S denotes the set of the acceptable solution, and this set is given by the constraint described by the radius α of the central composite plan of the experiments.
In the paper, the following form of the desirability function is used [11]: where d i is the desirability function related to the optimization's i-th criterion, w i denotes the weight of the i-th criterion, and m is the number of responses.Two types of the desirability function are applied: the smaller-the better (STB) and the nominal-the better (NTB) [11] as follows: STB : and NTB : where for the STB function, L denotes the acceptable target for the response ŷi , and U is an acceptable upper limit of the response, s is the exponent which sets the sharpness of the desirability function.For the NTB function, T denotes the target for the response ŷi , and L and U are acceptable lower and upper limits of the response, respectively.The responses ŷi in the presented formulation are obtained using the identified RSM models.Table 13 shows the parameters of the identified models' desirability functions.Table 13 shows that the most critical response is the first one-related to the deflection of the rack bar.
Using the first of Equation (3) with Equations ( 4) and ( 5) and Table 12, the desirability function was calculated for the responses given in Table 1 obtained during the experimental phase of the study.
The values of the desirability function for the plan of experiments results are given in Figure 9.

6.
HRCt3 Table 13 shows that the most critical response is the first one-related to the deflection of the rack bar.
Using the first of Equation (3) with Equations ( 4) and ( 5) and Table 12, the desirability function was calculated for the responses given in Table 1 obtained during the experimental phase of the study.
The values of the desirability function for the plan of experiments results are given in Figure 9.The classical optimization technique, namely the sequential quadratic programming (SQP) algorithm for constrained optimization problems, was applied in the first trial.One  The classical optimization technique, namely the sequential quadratic programming (SQP) algorithm for constrained optimization problems, was applied in the first trial.One hundred tests were carried out with a random selection of the starting point of the optimization procedure.The results are presented in Figure 10.It is visible that only for four out of one-hundred trials, the best value of the desirability function is non-zero with a maximum of 0.646.This poor result can be explained by approximating the gradient of the objective function using difference quotients.However, one must also consider that the desirability functions ( 4) and ( 5) are not differentiable in It is visible that only for four out of one-hundred trials, the best value of the desirability function is non-zero with a maximum of 0.646.This poor result can be explained by approximating the gradient of the objective function using difference quotients.However, one must also consider that the desirability functions (4) and ( 5) are not differentiable in the classical sense.In Table 14, the statistics for the SQP optimization process are presented.The results show that the process is not robust; the optimal design variables vector in uncoded values is x opt = [56.4% 1.6 mm 796 mm/min] T , which gives D opt = 0.646.The second best result is x = [56.4% 3.1 mm 793 mm/min] T which gives D opt = 0.557.The latter result is significant from the point of view of the industrial conditions of the process; additional considerations will also be presented in this work.Figure 11 shows slices of the desirability function hypersurface.The hypersurface of desirability has the form of narrow ribs with a wavy ridge, making it very difficult to solve the optimization problem using classical algorithms.In addition, for almost the entire set of admissible solutions, the values of the desirability function are equal to zero, which also justifies the inefficiency of the classic optimization algorithm SQP.Due to the presented features of the objective function, an evolutionary, floatingpoint coded genetic algorithm (GA) was used in the work.
A genetic algorithm (GA) [21] is an intelligent evolutionary procedure used to find the optimal solutions to problems based on natural genetics and natural selection principles.The foundation of this algorithm is the biologically inspired set of operations such as selection, combination or crossover, and mutation.In the present case, the individuals and the objective function are defined and coded at the real-coded strings.Then, the iterative process is carried up using the three operations, i.e., selection, crossover, and mutation.The selection operation is to choose individuals that start from a population generated randomly according to a probability distribution.The roulette wheel selection method is used to determine the selection of individuals.A crossover operation appears when two strings randomly picked from the populations exchange their substrings to create two new strings.The sum of individuals to apply the crossover operation is dominated by crossover probability, which is the ratio of all selected strings to the total population A genetic algorithm (GA) [21] is an intelligent evolutionary procedure used to find the optimal solutions to problems based on natural genetics and natural selection principles.The foundation of this algorithm is the biologically inspired set of operations such as selection, combination or crossover, and mutation.In the present case, the individuals and the objective function are defined and coded at the real-coded strings.Then, the iterative process is carried up using the three operations, i.e., selection, crossover, and mutation.The selection operation is to choose individuals that start from a population generated randomly according to a probability distribution.The roulette wheel selection method is used to determine the selection of individuals.A crossover operation appears when two strings randomly picked from the populations exchange their substrings to create two new strings.The sum of individuals to apply the crossover operation is dominated by crossover probability, which is the ratio of all selected strings to the total population size.The mutation operation assures the diversity of the population.It is an occasional random alternative in one or more string positions where the small random number mutes the value.After successive iterative generations, the population evolves gradually toward an optimal solution.Finally, the evolutionary process is completed till the termination criterion is reached.This study proposed RSM with the desirability function to define the relationship between parameters and responses.A practical method for combining RSM and an evolutionary strategy in the form of a genetic algorithm (GA) is that RSM is utilized to build a functional relationship between parameters and responses, and then GA is used to optimize the given fitness function composed by the desirability approach.
The following parameters of the applied GA were listed based on the study to get optimal solutions with the lowest computational effort: The constraint: the first generation is chosen randomly according to the constraint.
The optimization process is conducted 100 times, and the best results after each trial are collected in Figure 12.The optimization process is conducted 100 times, and the best results after each trial are collected in Figure 12.It is visible that for all 100 trials, the best value of the desirability function is non-zero, with an overall maximum of 0.646.In a few cases, local minima were found, or the optimization process was interrupted due to the lack of a significant change in the value of the objective function.Most often, the optimization process requires five iterations of the algorithm, which is a relatively small value, in particular for methods of computational intelligence, such as the genetic algorithm.In Table 15, the statistics for the evolutionary optimization process are presented.It is visible that for all 100 trials, the best value of the desirability function is nonzero, with an overall maximum of 0.646.In a few cases, local minima were found, or the optimization process was interrupted due to the lack of a significant change in the value of the objective function.Most often, the optimization process requires five iterations of the algorithm, which is a relatively small value, in particular for methods of computational intelligence, such as the genetic algorithm.In Table 15, the statistics for the evolutionary optimization process are presented.The results in Table 15 show that evolutionary optimization is a robust technique, especially for complex non-differentiable multidimensional fitness functions as the considered desirability function for the induction hardening process of the rack bar.
The global optimal design variables vector in uncoded values is x opt = [56.4% 1.6 mm 796 mm/min] T , which gives D opt = 0.646; this result agrees with the one obtained for the SQP algorithm considering the rounding used.One of the results for local minima is also close to the second-best result found with the SQP algorithm, namely x = [56.4% 3.0 mm 795 mm/min] T , which gives D opt = 0.556.However, it should be noted that the GA found more local minima, as can be seen in Figure 12, but similar results are shown to compare the results of the two algorithms.Figure 13 shows a histogram of the optimal values of the design variables obtained during 100 tests of the GA.It can be seen that the most significant variation in optimal values occurs for the third component of the vector of design variables, that is, the feed rate.For the first and second components of the vector of design variables, the optimal values cluster with less variation than for the third component.variables, that is, the feed rate.For the first and second components of the vector of design variables, the optimal values cluster with less variation than for the third component.

Confirmation Experiments and Discussion
The best and global optimal point obtained during the evolutionary optimization was found for the following values of the design variables: 56.4% of power, 1.6 mm of the distance of the hardening coil, and 804 mm/min of the feed rate.The second design variable was found close to its lower bound.Industrial conditions require a different approach to evaluating the presented optimal solution, which can be regarded with high probability as a global optimum.Preliminary tests have shown that the found optimal solution cannot be applied due to increased heating of the tool, which is a coil located close to the heated element, leading to its rapid wear, increasing the cost of the process.The original assumptions about the range of design variables included the possibility of fully controlling design variables without taking into account additional factors that arise in industrial conditions, especially in the conditions of mass production in an automatic cycle.Therefore,

Confirmation Experiments and Discussion
The best and global optimal point obtained during the evolutionary optimization was found for the following values of the design variables: 56.4% of power, 1.6 mm of the distance of the hardening coil, and 804 mm/min of the feed rate.The second design variable was found close to its lower bound.Industrial conditions require a different approach to evaluating the presented optimal solution, which can be regarded with high probability as a global optimum.Preliminary tests have shown that the found optimal solution cannot be applied due to increased heating of the tool, which is a coil located close to the heated element, leading to its rapid wear, increasing the cost of the process.The The result for the first test shows that the obtained values of the desirability function are not statistically equal to the calculated optimal value, which is mainly influenced by the results of the second and fourth experiments, for which the values of the desirability function are zero.The rest of the test results indicate that the process response requirements are statistically met at a significance level of 0.05.Moreover, process quality improvement is observed, mainly if the deflection response is analyzed as the most critical response.Uncertainties of the process mean that it is impossible to talk about finding the optimal parameters of the induction hardening operation, both in the sense of single-criteria and multi-criteria optimization.However, it could be concluded that optimizing the desirability function based on the RSM models increases the quality of the process.It can be noted that if the confirmation sample had a larger size, the process capability analysis could have helped in assessing its quality improvement.

Conclusions
This work addressed the constrained multi-response optimization of processing conditions in induction hardening to decrease mechanical deformation during the hardening operation and improve the process's predefined geometrical and mechanical properties.The RSM models were used to model the highly nonlinear relationships between inputs (power, distance, feed rate) and technological responses.The GA was applied to determine the quasi-optimal values of performances measured and technological inputs.
The following important conclusions can be drawn.The induction hardening process can be significantly improved.The quadratic models of responses of the process were significant according to the quality measures of the least-square approximation.The developed RSM models act effectively in the optimization process.The models proposed can be used in industrial applications to predict technological responses with acceptable accuracy.
The optimization problem, in which maximal deflection, depths of the hardening and hardness of teeth and the back side of the rack bar are objectives, and the desirability function is defined using them, is practical and realistic in the induction hardening process optimization compared to a single objective or simultaneously optimizing nine responses in the Pareto sense.However, other formulations are possible, such as one in which some criteria act as constraints; this may provide directions for further research.Industrial conditions, on the other hand, favor the use of desirability functions.Statistically significant lack of fit of some response functions indicates directions for further research related to using other empirical modeling tools, such as neural networks or kriging.However, this forces larger training data sets, which are difficult to obtain in industrial conditions.
The optimization algorithms recommended the following quasi-optimal combination of process parameters: 56.4% power, 3.1 mm distance of the hardening coil, and 793 mm/min of the feed rate when the industrial conditions were considered.Industrial conditions in the considered case forced the rejection of the possible global optimum of the desirability function, and the quasi-optimal solution was recommended, related to the local optimum, indicated approximately by both the classical optimization algorithm, such as SQP and the genetic algorithm.Uncertainties in empirical modeling also indicate that industrial conditions with limited experimental data (also with a lack of replication) make it difficult to determine optimal solutions.However, the confirmation experiments show that the quasi-optimal process is reliable and high-quality.The use of the nonparametric Wilcoxon test allowed the evaluation of a confirmatory experiment that provided a small sample of data, which is due to the industrial conditions of the research conducted.The test results confirm a statistically significant improvement in process quality in terms of the measures introduced.
The proposed hybrid approach in this paper can be considered an effective technique in some academic research and industrial applications to identify the optimal parameters in the induction hardening process and decrease production costs and time.
Industrial conditions should be carefully applied when solving optimization problems to meet multiple requirements.However, it can be concluded that the influence of

Figure 1 .
Figure 1.The automatic inductive hardening and tempering machine.

Figure 1 .
Figure 1.The automatic inductive hardening and tempering machine.

Figure 3 .
Figure 3.The general flow of the experimentation procedure.

Figure 4 .
Figure 4. Zones for the hardness measurements of the rack bar after induction hardening operati

Figure 3 .
Figure 3.The general flow of the experimentation procedure.

Figure 3 .
Figure 3.The general flow of the experimentation procedure.

Figure 4 .
Figure 4. Zones for the hardness measurements of the rack bar after induction hardening operation.

Figure 4 .
Figure 4. Zones for the hardness measurements of the rack bar after induction hardening operation.

Materials 2023 ,Figure 5 .
Figure 5.The setup of the hardness depth measurements.

Figure 5 .
Figure 5.The setup of the hardness depth measurements.

Figure 5 .
Figure 5.The setup of the hardness depth measurements.

Figure 6 .
Figure 6.Example of hardness profile from the tooth side.

Figure 6 .
Figure 6.Example of hardness profile from the tooth side.

Figure 7 .
Figure 7.The setup of the deformation measurements.
x s = [−0.1960.764 1.053] T is a saddle point because the eigenvalues of matrix B {−31.45, 36.93,269.55} are mixed in sign.The reduced model, according to Table

Figure 8 .
Figure 8. Contour plot for the slice x3 = 0.5 of the maximal deflection response hypers the full quadratic model; (b) the reduced model.

Figure 8 .
Figure 8. Contour plot for the slice x 3 = 0.5 of the maximal deflection response hypersurface y 1 : (a) the full quadratic model; (b) the reduced model.

Figure 9 .
Figure 9. Desirability values for the experimental data obtained using CCD.

Figure 9
Figure9shows that only three out of twenty experiments gave non-zero desirability values with a maximum of 0.421 for the seventh and 0.292 and 0.291 for the fourth and fourteenth experiments, respectively.The results indicate that there is room for improvement in the process.The classical optimization technique, namely the sequential quadratic programming (SQP) algorithm for constrained optimization problems, was applied in the first trial.One

Figure 9 .
Figure 9. Desirability values for the experimental data obtained using CCD.

Figure 9
Figure9shows that only three out of twenty experiments gave non-zero desirability values with a maximum of 0.421 for the seventh and 0.292 and 0.291 for the fourth and fourteenth experiments, respectively.The results indicate that there is room for improvement in the process.The classical optimization technique, namely the sequential quadratic programming (SQP) algorithm for constrained optimization problems, was applied in the first trial.One hundred tests were carried out with a random selection of the starting point of the optimization procedure.The results are presented in Figure10.

Materials 2023 ,
16,  x FOR PEER REVIEW 18 of 24 hundred tests were carried out with a random selection of the starting point of the optimization procedure.The results are presented in Figure10.

Figure 10 .
Figure 10.Best-found values of the desirability function for SQP optimization during N = 100 tests.

Figure 10 .
Figure 10.Best-found values of the desirability function for SQP optimization during N = 100 tests.

Figure 11 .
Figure 11.Slices of the desirability function hypersurface D(x 1 , x 2 , x 3 ): (a) x 3 = 0; (b) x 2 = 0; (c) x 1 = 0.Due to the presented features of the objective function, an evolutionary, floating-point coded genetic algorithm (GA) was used in the work.A genetic algorithm (GA)[21] is an intelligent evolutionary procedure used to find the optimal solutions to problems based on natural genetics and natural selection principles.The foundation of this algorithm is the biologically inspired set of operations such as selection, combination or crossover, and mutation.In the present case, the individuals and the objective function are defined and coded at the real-coded strings.Then, the iterative process is carried up using the three operations, i.e., selection, crossover, and mutation.The selection operation is to choose individuals that start from a population generated

Figure 12 .
Figure 12.Best-found values of the desirability function for GA optimization during N = 100 tests.

Figure 12 .
Figure 12.Best-found values of the desirability function for GA optimization during N = 100 tests.

Figure 13 .
Figure 13.Histograms of optimal design variables for GA optimization results.

Figure 13 .
Figure 13.Histograms of optimal design variables for GA optimization results.

Table 1 .
The CCD plan of experiments with process responses.

Table 1 .
The CCD plan of experiments with process responses.

Table 2 .
Analysis of variance for the mathematical model of maximal deflection response y 1 .

Table 3 .
Analysis of variance for the mathematical model of the hardening depth of teeth y 2 .

Table 4 .
Analysis of variance for the mathematical model of the hardening depth on the back side y 3 .

Table 5 .
Analysis of variance for the mathematical model of the hardness of the teeth in the first zone y 4 .

Table 6 .
Analysis of variance for the mathematical model of the hardness of the teeth in the second zone y 5 .

Table 7 .
Analysis of variance for the mathematical model of the hardness of the teeth in the third zone y 6 .

Table 8 .
Analysis of variance for the mathematical model of the hardness of the back-side in the first zone y 7 .

Table 9 .
Analysis of variance for the mathematical model of the hardness of the back-side in the second zone y 8 .

Table 10 .
Analysis of variance for the mathematical model of the hardness of the back-side in the third zone y 9 .

Table 11 .
Statistical characteristics of regression models for responses of the process.

Table 12 .
The reduced models of responses.

Table 13 .
Parameters of the desirability functions for the identified response models.

Table 14 .
Statistics of the SQP optimization results for the desirability function values.

Table 15 .
Statistics of the evolutionary optimization results for the desirability function values.