1 Introduction

Computer numerical controlled (CNC) milling is an essential metal cutting technique among various machining processes in the modern era of manufacturing. CNC milling not only makes the milling process fully automated but also enhances machining time, reduces process variations, improves product quality, and enriches the overall productivity of the manufacturing companies. For CNC milling, end mill is the most exploited among various milling cutters due to its ability of high-speed cutting of metal with minimum surface roughness in a single pass [1]. CNC end milling is being used significantly in different manufacturing sectors, such as aerospace, automotive, electronics, jewelry, and bioinstrumentation industries. CNC end milling is used for making different geometrical shapes and holes in a metallic workpiece during milling, profiling, contouring, slotting, counter-boring, drilling, reaming applications, etc. [2]. Aluminum alloy is the mostly explored material for end milling, which has more than 90% pure aluminum. It has high strength and ductility, corrosion resistance, weldability, machinability, and formability. It is used as an important material for vehicle bodies, refrigerated trucks, cold storage rooms, anti-skid flooring, manufacturing of mobile homes, residential siding, sheet metal work, rain carrying goods, etc. [3]. With optimum settings of end milling parameters, it is possible to achieve good surface quality and higher metal removal rate (MRR) for the aluminum alloy. For that matter, tool diameter (TD), spindle speed (SS), feed rate (FR), and depth of cut (DOC) could be the most important process variables. The surface roughness (Ra) is the primary machining attribute since most of the manufacturing companies try to maintain better surface quality for the machined parts. Therefore, Ra determines the manufacturing cost and quality of the engineered products [4]. Surface texture, fatigue resistance, and heat transmission of manufactured products are greatly influenced by Ra. Surface quality is also influenced by the abovementioned machining parameters of end milling. On the contrary, MRR is determined by the volume of removed metal and required machining time. MRR could affect the cost of manufacturing largely. When the combined effect of MRR and Ra is studied, the cutting process optimization becomes more complicated [5]. The solidity and life of the cutting tools could be greatly influenced by another process response known as cutting forces for the end milling process. Machining errors could be seen if the cutting forces are not considered while optimizing the process [6] [7]. Therefore, it is obligatory to study the effect of various process parameters on Ra, MRR, and cutting forces collectively for the CNC end milling process. The objective of this research is to determine the ideal parameter settings, which could yield high MRR and low Ra and cutting forces for the end milling of aluminum alloy (3105). This type of multi-response process optimization has never been accomplished on the Proxxon CNC milling machine, which is the focus area of this work.

To accomplish the above task successfully, the primary analysis is performed in two stages. In the first stage, Taguchi’s orthogonal array (OA) L16 (4 × 4) is used for designing the experimental study for end milling of AA3105 alloy. Kohonen’s self-organizing map (KSOM) is applied, which reveals the correlations among the process responses. Thereafter, grey relational analysis (GRA) and VIseKriterijumska Optimizacija I Kompromisno Resenje (VIKOR) methods are employed as multi-response optimization techniques, which can simultaneously optimize MRR, Ra, and cutting forces by obtaining single grey relational grade (GRG) and VIKOR index score respectively. Furthermore, the analysis of variance (ANOVA) is conducted to understand the most influential process parameter and its percentage contribution. In the second stage, the Bayesian regularized neural network (BRNN)-assisted beetle antennae search (BAS) algorithm is developed and used as surrogate-assisted optimization technique for the end milling of AA3105. This sequential BRNN-BAS model accurately estimates the process parameter settings for maximum MRR, and minimum Ra and cutting forces. The BRNN is used as a surrogate module and BAS is used as the optimization module in this study. Another variant of BAS algorithm is introduced to elaborate the analysis, which exploits multiple regression model as the surrogate. The rest of this article is arranged in the following order, in the “Related work” section, the related literature are cited, which covers the papers of end milling of alloys and composites, deep learning-based predictive modeling of manufacturing processes, and surrogate-assisted data-driven optimization. The “Experimental setup and materials” section demonstrates the details of material, experiments on Proxxon milling machine, and correlational analysis among input and output variables. BRNN-BAS-based data-driven surrogate-assisted optimization technique is presented in the “Research methodologies” section along with the GRA, VIKOR method, and regression-BAS algorithm. Computational studies, analysis of results, and discussions are presented in the “Results and discussions” section. Finally, the “Conclusions” section concludes this work and shows the future scopes.

2 Related work

2.1 CNC milling process optimization

Due to the tremendous advancements in material research and increasing complicacies in part geometries, milling cutters are evolving in faster rate. End mill tools are being applied heavily in various high-speed metal removal processes. Recently, different types of milling tools are being investigated by manufacturing researchers, which are flat/square end mill, ball end mill, fillet mill, chamfer mill, face mill, etc. Li and Zhu [8] utilized a flat end mill tool for cutting force prediction modeling on an improved automatic five-axis machine, which is further coupled with the cutting effects. Chao and Altintas [9] used the mechanics and dynamics of the ball end milling process on a 5-axis milling machine. Zhu [10] tried the fillet mill for the feature-based modeling with parametric design of the cutting process. Gao et al. [11] utilized various chamfered mills for slot milling of aluminum alloy 7075, where the tool wear and surface roughness are considered as the performance characteristics. Zhenyu et al. [12] studied the face milling with an algorithm and considered the different factors of Ra for estimation. In this paper, the square end mill is used as the cutter as displayed in Fig. 1. This is the mostly common end mill tool available for the CNC machine. Optimal parametric design for CNC milling has always been an area of interest to the researchers since the past two decades [13] [14] [15] [16] [17]. During this time, the methodologies for the multi-response milling processes have been advanced considerably. Mukherjee and Ray [18] presented a comprehensive study on the methodologies developed for the multi-response CNC milling. Yusup et al. [19] portrayed a detailed survey on the evolutionary algorithms (EA) and bio-inspired techniques for parametric designs of machining processes during 2007–2011. The authors concluded that the genetic algorithm (GA)-based optimization approaches are heavily exploited and the Ra is the most explored performance criterion. Furthermore, Chandrasekaran et al. [20] reviewed 142 articles based on the soft-computing and machine learning techniques for multi-response CNC machining. The authors stated that despite using the artificial neural network (ANN) and EA-based techniques, few more improvements are expected as future goal depends upon the data acquisition cost, the data filtration for noise, and the feature extraction from data. Kadirgama et al. [21] has proposed an ant colony optimization (ACO)-based technique to obtain optimum surface roughness for milling operation.

Fig. 1
figure 1

Six-mm flat/square end mill tool

Full size image

The technique is developed using the response surface method (RSM) and ACO. The most crucial variables are determined as cutting speed, feed rate, axial depth, and radial depth. Yildiz [22] proposed a hybrid algorithm to obtain optimal milling parameters. The results were successfully compared with GA, PSO, ACO, and immune algorithms. Arnaiz-González et al. [23] demonstrated the ball end milling process models using the multi-layer perceptron (MLP) and radial basis function (RBF). The RBF is shown to obtain better predictive model than the MLP with higher precision. Xiang and Zhang [24] depicted a prediction model based on the back propagation (BP) network and support vector machine (SVM) for milling process modeling, and proposed an optimization technique using the SVM and NSGA-II. Sarıkaya et al. [25] used the vibration signals, cutting force, and surface roughness to model the multi-response process of CNC milling of AISI 1050 steel. Das et al. [26] studied the grey fuzzy technique for the optimal level of process variables for CNC milling of Al–4.5% Cu–TiC metal matrix composites (MMCs). The cutting force (Fc) and surface roughness (Ra and Rz) are considered as performance characteristics. Zhou et al. [27] proposed an integrated multi-response technique based on GRA, RBF, and particle swarm optimization (PSO) for the ball end milling. While compared with the classical GRA, it produces continuous space for optimality. The authors considered some design variables such as inclination angle, cutting speed, feed, surface roughness, and compressive residuals stress. Khorasani and Yazdi [28] proposed the Ra monitoring system for the milling process considering cutting speed, feed per tooth, cut depth, type of materials, and coolant fluid as the decision variables and mechanical vibrations, white noise, and Ra as the process responses. Thereafter, the testing and verification procedures are utilized to achieve higher accuracy. Das et al. [29] performed the CNC milling of alumina green ceramic compact, where spindle speed, XY speed, Z speed, and depth of cut are considered as the process parameters and surface roughness is considered as the performance indicator. The authors developed a regression model from the experimental results and performed GA-based optimization of the Ra. Gaikhe et al. [30] experimented with the cutting forces in helical ball end milling of Inconel 718 and developed a GA by using cutting force regression model to obtain optimum cutting parameters. The proposed work achieved stability of machines, reduced the power consumption, and minimized the cost for machining. Kaushik et al. [31] has developed a statistical model to predict the temperature rise as a function of helix angle, the radial rake angle of the cutting tool. Cutting speed, feed rate, and axial depth of cut are considered as the process variables. RSM is used for optimization of high-speed end milling of aluminum Al 7068. Recently, Li et al. [17] proposed a multi-objective model for energy-efficient CNC milling with minimum production time and portrayed depth of cut and width of cut as the most influential parameters.

As of now, the following issues could be pointed out from the above citations,

  • Multi-response models are being considered for the CNC milling, where Ra, MRR, cutting forces, vibration, energy consumptions, etc., are considered as the process responses [20].

  • Some multi-response models are solved using the GRA or fuzzy-GRA-based single-objective techniques. The multiple regression, RSM, MLP, RBF, etc., are used for predictive modeling of process response [32].

  • Metaheuristic algorithms such as the GA, PSO, and ACO are used as the optimization tools. Mostly mathematical formulations are being used as the objective functions. Limited studies could be seen, where response surface or regression equations are used as the objective functions [21].

  • In true sense, the surrogate-assisted optimization (the deep learning and metaheuristic algorithms in combined form) has not yet been utilized for the manufacturing process optimization.

2.2 Surrogate-assisted optimization

Lately, the data-driven predictive or optimization models are being exploited more by the researchers for the CNC milling process. However, the combination of predictive and optimization models has not yet been studied completely, which is also known as data-driven surrogate-assisted optimization technique. It is an important state-of-the-art area for the machine learning and optimization researchers [33]. When the metaheuristic techniques are being used for the multi-response manufacturing process modeling, the best practice is to utilize the data-driven surrogate models as the fitness functions. This approach eliminates the requirement of the complex mathematical function, expensive empirical or numerical modeling for fitness evaluation. It often facilitates the use of traditional or existing optimization algorithms. Surrogate-assisted optimization could also be classified as black-box optimization when there would be little to no information available about the process or problem considered [34]. Yet this approach is highly capable of approximating the functional relationships among process variables based on the sampled data obtained using design of experiment (DOE) techniques [35]. The exactness of the surrogate training could be an important issue for the data-driven surrogate models if global optimality is not guaranteed. Therefore, the mean square error (MSE), root mean square error (RMSE), or mean absolute error (MAE) could be used as the performance metrics. The lower is the error score, the better is the accuracy of the model. Once the surrogate model is trained, an appropriate metaheuristic algorithm, e.g., GA, PSO, and ACO, could be employed to find optimal set of the process variables [36]. Data-driven surrogate-assisted optimizers are substantially prompt and efficient; therefore, these are computationally inexpensive. DOE-based tools, such as the central composite design (CCD), Box-Behnken design (BBD), D-optimal design (DOD), Latin hypercube sampling (LHS), full factorial design (FFD), and orthogonal array design (OAD), are generally used to define the experimental or trial sample points as initial solution(s). These are also used to define the input data to the surrogate models. DOE methods generally maximize the process information obtained from the restricted number of trial runs [37]. As mentioned earlier, some of the surrogate models already exploited for the machining process modeling are RSM, Gaussian process (GP), decision tree, RBF, SVM, etc. [32]. Recently, surrogate-assisted optimization techniques are being practiced in various engineering and technological research. Sun et al. [38] developed a two-stage surrogate-assisted PSO algorithm by incorporating a global and some local surrogate functions. The proposed technique is tested with some popular unimodal and multimodal problems from literature. In another research [39], the authors demonstrated the online- and offline-based classifications of the surrogate-assisted optimization techniques and developed an EA to optimize the offline data-driven trauma system. The CPU time is reduced using a multi-fidelity surrogate management strategy. Haftka et al. [40] performed a survey based on the surrogate-based global optimization with a focus on the balance between the exploration and exploitation search and kriging-based surrogate models are reviewed primarily. Furthermore, Sun et al. [41] studied the combined effect of the surrogate-assisted PSO and a surrogate-assisted social learning (SL-PSO) algorithm. SL-PSO worked on exploration and PSO worked on exploitation. The proposed method is successfully tested on the benchmark problems. Allmendinger et al. [42] presented another survey and discussed the complexities in the surrogate-assisted multi-objective optimization. The authors found these complexities from the different real-world problems and analyzed. This study pointed out multiple future scopes and demonstrated the applicability of the surrogate-assisted optimization in industrial settings. Chugh et al. [43] proposed a kriging surrogate-assisted reference vector-guided EA (RVEA) and tested on some benchmark problems. Jin et al. [44] considered five real-world cases of blast furnace optimization, trauma system design optimization, optimization of fused magnesium furnaces, optimization of airfoil design, and air intake ventilation system optimization where surrogate-assisted optimizers have shown promising solutions.

Following points could be extracted further,

  • Surrogate-assisted optimization is useful for computationally expensive and many-objective optimization problems. These problems are based on real-world data, which are not readily available in the literature. Therefore, this area of research is less explored [44].

  • Different surrogate models are available in literature such as GP, kriging, and RBF for the surrogate-assisted optimization. Many other predictive models, such as MLP, BRNN, SVM, decision tree, Gaussian kernel, and deep belief network, could also be explored for the surrogate-assisted optimization [42].

  • Metaheuristic techniques are being evolved every day. Recently, researchers exhibit great enthusiasm in transforming bionic phenomenon into the computer algorithm. As a result, various new and promising methods have appeared, such as the grey wolf optimization (GWO), African buffalo optimization (ABO), beetle antennae search (BAS) algorithm, ant-lion optimization (ALO), Harris hawks optimization algorithm (HHO), the grasshopper optimization algorithm (GOA), and the multi-verse optimization algorithm (MVO) [45,46,47,48,49]. These latest methodologies are capable enough to outperform popular EAs and these are not much exploited for the surrogate-assisted optimization. Hence, a missing link exists between the deep learning and bio-inspired metaheuristic research, which is an innovation issue. This further shows immense research scope on the hybridization of methodologies.

  • Many different real-world problems are considered for the surrogate-assisted optimization. However, manufacturing or machining processes are not very popular yet for the said purpose. The reason could be that the manufacturing researchers have not yet shown their interest in the said approaches, or some communication gap exists among machine learning, optimization, and manufacturing researchers.

The aim of this article is to address these abovementioned issues. For that matter, a Bayesian regularized neural network (BRNN) is used for the surrogate modeling and the beetle antennae search (BAS) algorithm is used as the optimization technique. Both the techniques are novel and never employed in manufacturing process optimization. Hence, the proposed BRNN-BAS is verified further with the end milling cutting of aluminum alloy.

3 Experimental setup and materials

The experiments are carried out on Proxxon FF 500/BL 3-Axes CNC milling machine, manufactured by Proxxon, Germany. It has brushless direct drive with 50 μm precision, which is suitable for academic projects or small businesses. It has double roller bearing recirculating ball spindles at all 3-axes. The spindle speed varies in the range of 200–4000 rpm. It has large traverse area (X 290 mm, Y 100 mm, and Z 200 mm). In this machine, milling head could be pivoted to the left and right by 90°, with running condenser motor. The experimental setup is depicted in Fig. 2a. Four milling tools are considered (Fig. 2b) based on the spiral design according to DIN 844 and made of the high-speed steel (HSS-Co5) 5% cobalt. The cutting force data acquisition (DAQ) system consists of a KISTLER dynamometer 9257B (6-DOF), a KISTLER multichannel charge amplifier 5017B, a KISTLER DAQ system 5697A1, and the DynoWare software. To measure the surface quality of the machined workpieces, a ZEISS Handysurf E-35B is used.

Fig. 2
figure 2

a Proxxon CNC milling hub and b various end milling tools

Full size image

3.1 Materials

AA3105 alloy is used as the test workpiece (90 × 140 × 20 mm3) for the end milling cutting (Fig. 3). AA3105 is primarily used in sheet metal work and manufacture of mobile homes, residential siding, and rain carrying goods in sub-zero temperature. It is perfectly suitable for the climate of Nordic Europe. AA3105 portrays good machinability property. Percentages of weight in chemical composition of the AA3105 are Al 98.56%, Mn 0.716%, Fe 0.38%, Zn 0.128%, Cu 0.118%, Cr 0.081%, and Pb 0.006%. Table 1 shows the mechanical properties of the alloy. The workpieces are positioned and clamped with fixtures supplied by Proxxon and KISTLER. Many independent process variables are involved in this experiment. To determine the most influential factors, rigorous experimentations were carried out on AA3105 alloy. Specific ranges for the process variables are determined from machine handbook and empirical studies. Cutting force data sample is portrayed in Fig. 4, which is recorded in an interval of 0.01 s for a 10-s cycle.

Fig. 3
figure 3

Sample AA3105 workpiece

Full size image
Table 1 Mechanical property of AA3105
Full size table
Fig. 4
figure 4

Sample of cutting force data

Full size image

3.2 Experiments

TD, SS, FR, and DOC are selected as most influential process parameters for the experiments carried out. For Taguchi’s OAD, the level of each factor is opted based on the operating range. Parameters are set with four levels and portrayed in Table 2.

Table 2 Milling parameters and their levels for L16 orthogonal array (OA)
Full size table

In this study, MRR, Ra, and cutting forces in three directions (Fx, Fy, and Fz) are considered as process responses. Taguchi’s L16 OAD is considered, which consists of sixteen experimental runs. These experiments are conducted on fresh workpieces for three repeated runs each for every setting and the average values are recorded. The experimental results are depicted in Table 3. The main effect plots for the responses are portrayed in Fig. 5, which points out the following observations,

  • The MRR is directly correlated and proportionate with the TD and DOC. This means the higher is the TD and DOC scores, the better is the MRR score, whereas the SS and FR are inversely related to the MRR, i.e., the MRR score is higher at the low-level settings of the SS and FR. MRR is dependent mostly on DOC and FR. Hence, a higher DOC indicates lower FR to achieve a level of MRR. Furthermore, the low FR demands lower SS in order to improve the tool life [15]. Therefore, the ideal setting for the MRR is TD4-SS1-FR3-DOC4 which produces better mean value for MRR (Fig. 5).

  • For the high value of mean Ra, high scores of the TD, DOC, and SS are required and lower value of DOC is expected. For better-quality surface, lower DOC and higher SS are desirable in CNC milling [2]. The desired setting for Ra is TD4-SS4-FR2-DOC4 (Fig. 5).

  • Three cutting force components (Fx, Fy, and Fz) change their values abruptly with small changes on parameters. For the Fx, TD4-SS4-FR2-DOC1 would be the ideal setting, and for the Fy, TD2-SS2-FR1-DOC1 yields best result, whereas for the Fz, TD4-SS2-FR1-DOC1 is the optimal combination. It could be concluded that determining the exact levels of the process parameters is a combinatorial task while attaining the near-optimal force components together with the MRR and Ra.

Table 3 Experimental results for MRR, Ra, Fx, Fy, and Fz in Taguchi’s L16 OAD
Full size table
Fig. 5
figure 5

Individual main effect plots for MRR, Ra, Fx, Fy, and Fz

Full size image

3.3 KSOM analysis

The Kohonen’s self-organizing map (KSOM) is an unsupervised deep learning tool, which presents the high-dimensional (n-D) data using a two-dimensional (2D) map [50]. This presentation saves the topological information of the original data in the new data. This graphical 2D presentation enables the users to visualize data correlations [51] [52]. KSOM trains the network using an input data point x. Using competitive learning, the winning neuron c with the weight vector wc is selected using,

$$ \left|\left|x-{w}_c\right|\right|={\min}_i\left(\left|\left|x-{w}_i\right|\right|\right) $$
(1)

The weight vector wc is updated to match the data point x. Thereafter, the weights of the neurons in the close neighborhood of c are also updated, i.e., the nearby neurons are moved towards c. The update expression is demonstrated as,

$$ {w}_i\left(t+1\right)={w}_i(t)+{h}_{ci}(t)\left[x(t)-{w}_i(t)\right] $$
(2)

where t is the count of iterations, x(t) is the input data point x during t, and hci(t) is the neighborhood function (decreasing type Gaussian function), defined as,

$$ {h}_{ci}(t)=\alpha (t)\times {e}^{\left(-\frac{\left|{r}_c-{r}_i\right|\mid }{2\sigma {(t)}^2}\right)} $$
(3)

where α(t) is the learning rate. KSOM algorithm reduces the dimensionality of the data to present it in 2D and reveals the hidden pattern in data.

Figure 5 previously demonstrated the relationships among inputs and responses. Hence, the KSOM is employed to analyze the correlations among all the five responses. In order to obtain high quality of KSOM topology structure, the quantization error and the topological errors are kept at the minimum level (≃ 0). The size of the map is considered as 19 × 17. The KSOM visualization of data is presented in Fig. 6. The KSOM maps are the distribution of colors (from lighter shed to darker shed) with the actual values of variables. It could be stated that the MRR and Ra are inversely related as the MRR map shows light color at the top-right corner of the map and the Ra depicts dark color distribution at the same area of the map. This phenomenon can be scientifically explained. When MRR is high, a higher number of atomic layers are chipped off the surface of the material. Due to the higher speed of the chipped particles, the atoms of the material surface are chipped off evenly in large chunk. This leads to the surface of higher quality [28]. Force components along the x-axis and y-axis show very similar nature of color patterns. However, third force component along z-axis shows different pattern of colors. It is also notable that the MRR and Ra obtain low scores when forces are at higher side and vice versa. This study implies that the responses are conflicting in nature, and to optimize the multi-response end milling process, concurrent optimization of all the responses is required to obtain the near-optimal levels of the end milling process parameters.

Fig. 6
figure 6

Visualizations of correlations among inputs and responses for end milling of AA3105 alloy

Full size image

4 Research methodologies

The proposed multi-response optimization is performed using four different approaches, such as GRA, VIKOR methods, BRNN-assisted BAS, and regression-assisted BAS algorithms. Furthermore, the comparative analyses are presented based on these methods, and the best solution is picked up. In this section, the methodologies are discussed.

4.1 Taguchi’s OAD coupled with grey relational analysis

Taguchi’s OAD is suitable for single-response design. For multi-objective problems, the GRA is suitable [53], which has the ability to exploit Taguchi’s OAD and estimate the process responses using single grey relational grade (GRG). Steps of GRA are as follows:

  1. Step 1:

    The data are normalized to reduce inconsistency and restricted in the range {0, 1}. When the performance objective is to be minimized, non-beneficial (Eq. 4) rule is applied; else, beneficial (Eq. 5) rule is applied,

$$ {y}_i^{\ast }(x)=\frac{y_i^0{(x)}_{\mathrm{max}}-{y}_i^0(x)}{y_i^0{(x)}_{\mathrm{max}}-{y}_i^0{(x)}_{\mathrm{min}}} $$
(4)
$$ {y}_i^{\ast }(x)=\frac{y_i^0(x)-{y}_i^0{(x)}_{\mathrm{min}}}{y_i^0{(x)}_{\mathrm{max}}-{y}_i^0{(x)}_{\mathrm{min}}} $$
(5)

where, i∈[1, m] and x∈[1, N], m is the number of experimental runs, and N is the number of responses. yi0(x)max and yi0(x)min are the largest and smallest values of yi0(x) is the original data.

  1. Step 2:

    Compute grey relational coefficient (GRC) using Eq. (6),

$$ {\varepsilon}_i(x)=\frac{\delta_{\mathrm{min}}-\varepsilon \times {\delta}_{\mathrm{max}}}{\delta_i^0(x)-\varepsilon \times {\delta}_{\mathrm{max}}} $$
(6)

where δi0 (x) = yi0 (x)-yi*(x), the deviation coefficient.

  1. Step 3:

    Calculate grey relational grade (GRG) using Eq. (7),

$$ {\gamma}_i=\frac{1}{N}\times \sum \limits_{k=1}^N{\varepsilon}_i(x)\kern0.5em $$
(7)

GRG depicts the overall quality index (transformed single response). The values of GRG determine the ranking of experimental runs and obtain near-optimal set of variables.

  1. Step 4:

    Calculate the analysis of variance (ANOVA) to find out the sensitivity of the variables to the design process at 95% confidence level and obtain the response table and main effect plot. This further depicts ranks based on delta statistics, which compare the relative magnitude of effects. The delta statistic shows the difference between the largest and the smallest average (max-min) for each variable. It finally indicates the most sensitive variables to the design process.

  2. Step 5:

    Once the optimal levels of process parameters are selected, the improvement of GRG is measured using the confirmatory test. The predicted GRG could be computed using,

$$ \hat{\gamma}={\gamma}_m+\sum \limits_{i=0}^n\left({\gamma}_i-{\gamma}_m\right) $$
(8)

where γm is the sum of mean GRG, γi is the mean GRG at the near-optimal level of the individual process parameter, and n is the number of input variables.

4.2 Taguchi’s OAD coupled with VIKOR method

The VIKOR method is a multi-criteria decision-making (MCDM) technique, which was originally developed in reference [54] to solve the decision problems with conflicting objectives. It is based on a decision matrix with alternatives (rows) and responses (columns). The conflicting correlations are investigated, and solutions are obtained, which are closest to the ideal. Furthermore, the alternatives (rows) are evaluated based on the verified conditions. The VIKOR method assigns ranks to the alternatives based on the distance to the ideal. The VIKOR method could be applied on Taguchi’s OAD to transform the multi-response process into a single-response process based on VIKOR index. The steps are as follows,

  1. Step 1.

    An initial decision matrix is created: where Ai represents ith alternative, i=1, 2, 3, … 푚; Cxj represents the jth criterion, j = 1, 2,.........,n, and xij is the individual performance of an alternative.

$$ D=\left[\begin{array}{ccc}{C}_{x_1}& {C}_{x_2}& .\kern0.5em .\kern0.5em {C}_{x_n}\\ {}{x}_{11}& {x}_{12}& .\kern0.5em .\kern0.5em {x}_{1n}\\ {}\begin{array}{c}{x}_{21}\\ {}.\\ {}{\mathrm{x}}_{m1}\end{array}& \begin{array}{c}{x}_{22}\\ {}.\\ {}{x}_{m2}\end{array}& \begin{array}{c}\begin{array}{ccc}.& .& {x}_{2n}\end{array}\\ {}\begin{array}{ccc}.& .\kern1.25em & .\end{array}\\ {}\begin{array}{ccc}.& .& {x}_{mn}\end{array}\end{array}\end{array}\right] $$
(9)
  1. Step 2.

    The normalized decision matrix can be expressed as follows: F = [fij]m × n. Here, \( {f}_{ij}=\frac{x_{ij}}{\sqrt{\sum_{i=1}^m{x}_{ij}^2}} \), i = 1, 2, 3, …m; xij is the performance of Ai alternative with respect to the jth criterion.

  2. Step 3.

    For maximization problem, the ideal solution \( {f}_j^{\ast } \) and the negative ideal solution \( {f}_j^{-} \) are selected as \( {f}_j^{\ast }={\max}_i{f}_{ij} \), \( {f}_j^{-}={\min}_i{f}_{ij} \), and for minimization problem, the ideal solution \( {f}_j^{\ast } \) and the negative ideal solution \( {f}_j^{-} \) are selected as \( {f}_j^{\ast }={\min}_i{f}_{ij} \), \( {f}_j^{-}={\max}_i{f}_{ij} \)

  3. Step 4.

    Calculation of utility measure (Si) and regret measure (Ri): The utility measure and the regret measure for each alternative are given in (Eqs. 10–11),

$$ {S}_{ij}={w}_j\frac{\left({f}_j^{\ast }-{f}_{ij}\right)}{\left({f}_j^{\ast }-{f}_j^{-}\right)},{S}_i={\sum}_{j=1}^n{w}_j\frac{\left({f}_j^{\ast }-{f}_{ij}\right)}{\left({f}_j^{\ast }-{f}_j^{-}\right)}\kern0.5em $$
(10)
$$ {R}_i=\underset{j}{\min}\left({S}_{ij}\right)=\underset{j}{\min}\left(\ {w}_j\frac{\left({f}_j^{\ast }-{f}_{ij}\right)}{\left({f}_j^{\ast }-{f}_j^{-}\right)}\right) $$
(11)

where Si, and Ri represent the utility measure and the regret measure, and wj is the weight of the jth criterion, expressing the relative importance of each criterion. wj is selected using expert opinion or other MCDM techniques.

  1. Step 5.

    Computation of VIKOR index Qi (ith alternative VIKOR value, i = 1, 2, 3 … m) for the ith alternative is done using,

$$ {Q}_i=\frac{v\left({S}_i-{S}^{-}\right)}{\left({S}^{\ast }-{S}^{-}\right)}+\frac{\left(1-v\right)\left({R}_i-{R}^{\ast}\right)}{\left({R}^{\ast }-{R}^{-}\right)} $$
(12)
$$ {S}^{\ast }={\max}_i{S}_i,{S}^{-}={\min}_i{S}_i $$
(13)
$$ {R}^{\ast }={\max}_i{R}_i,{R}^{-}={\min}_i{R}_i $$
(14)

The term v is introduced as the weight of the maximum group utility. It ranges between [0, 1]. Ideally, it is set to 0.5. The VIKOR index is ranked based on the ascending order to find the best alternatives. Furthermore, the ANOVA is performed to find out the sensitivity of the variables to the design process at 95% confidence level and obtain the response table and main effect plot. Finally, delta statistics are obtained, and the most influential parameters are detected.

4.3 BRNN-BSA technique

The data-driven surrogate-assisted optimization technique is discussed here, which approximates the optimal combination of input variables for the end milling process. This technique can estimate the process responses accurately. The BRNN is opted for the predictive modeling and the BAS algorithm is considered as the optimization technique. BRNN-BAS is a sequential method, where BRNN executes first. Once BRNN training is completed, the best BRNN model is obtained with minimum RMSE score. Thereafter, the BAS optimization module executes.

4.3.1 BRNN model

BRNN is conceptualized on Bayes’ rule, which considers the relationship between the probability of any process and prior knowledge of the process. The output responses could be generated using this relationship, which is posterior. Foresee and Hagan [55] expressed the probability density function of weights as,

$$ P\left(w|d,\alpha, \beta, M\right)=\frac{P\left(D|w,\beta, M\right)P\left(w|\alpha, M\right)}{P\left(D|\alpha, \beta, M\right)} $$
(15)

where D is the data, M is the type of ANN (cascade forward multi-layer perceptron), w is the vector of network weights. P(w|α, M) denotes the prior knowledge of the weights. P(D|w, β, M) is the likelihood and P(D|α, β, M) is normalization factor. When noise in D and P(w|α, M) are Gaussian, then the likelihood and normalization factor are expressed as,

$$ P\left(D|w,\beta, M\right)=\frac{1}{Z_D\left(\beta \right)}{e}^{-\beta {E}_D} $$
(16)
$$ P\left(w|\alpha, M\right)=\frac{1}{{\mathrm{Z}}_w\left(\alpha \right)}{e}^{-\alpha {E}_w} $$
(17)

where ZD(β) = (π/β)n/2 and Zw(βα) = (π/α)N/2. ED is the sum of squared errors for data and EW is the sum of squares for weights. Hence, Eq. (15) becomes,

$$ P\left(w|d,\alpha, \beta, M\right)=\frac{\frac{1}{Z_D\left(\beta \right){Z}_w\left(\alpha \right)}{e}^{-\left(\beta {E}_D+\alpha {E}_w\right)}}{P\left(D|\alpha, \beta, M\right)} $$
(18)
$$ =\frac{1}{Z_F\left(\alpha, \beta \right)}{e}^{-\left(\beta {E}_D+\alpha {E}_w\right)} $$
(19)

Optimal weights maximize the posterior probability density function P(w|D, α, β, M), which is equivalent to minimization of regularized objective function f = (β.ED + α.Ew). This objective function restricts the weights and biases to be small enough, which produces smooth responses by eliminating over fitting. α and β are governed by Bayesian regularization. BRNN is a probabilistic ANN model. The input variables are generally treated as probability density functions to the hidden layer [56]. This method helps eliminating the uncertainty regarding the network weights. Therefore, the output obtained is spanned over a new distribution, which improves the prediction accuracy. The BRNN considered in this study has a multi-layer perceptron (MLP) network structure with input layer, one hidden layer, and output layer. Figure 7 shows the BRNN schematic. Various probability density functions are shown for the input and output.

Fig. 7
figure 7

BRNN schematic diagram

Full size image

4.3.2 BAS algorithm

Beetle antennae search (BAS) algorithm is a recently proposed technique, which is inspired from the odor-sensing mechanism of the beetles using their long antennae (Fig. 8) [57]. The antennae work as sensors with complex mechanism. Fundamental functions of such sensors are to follow the smell of food, or to sense the pheromone produced by the potential opposite gender for reproduction. The beetle moves its antennae in a particular direction to sense the smell of food or mates. This movement is random in the neighborhood area and directed according to the concentration of smell. This implies that the beetle turns to the right or the left depending upon the higher concentration of smell or odor data gathered by the antennae sensors.

Fig. 8
figure 8

Beetle search procedure based on odor-sensing mechanism using antennae

Full size image

Wang and Chen [58] described the BAS algorithm as,

Algorithm 1:

BAS algorithm for global minimum searching

figure a

The BAS algorithm is single-solution based, which is similar to the simulated annealing (SA) algorithm. The BAS starts with a randomly generated beetle with position vector xt at tth time (t = 1, 2, …, n) and the position is evaluated using the fitness function f(x), which determines the smell or odor concentration [59]. This fitness function could be a mathematical objective function, i.e., min f(x), or it could be an approximation model. For the purpose of initial solution generation, Latin hypercube sampling (LHS) is preferred, which is a popular method to sample the random solutions [60] [61]. In this paper, a modified form of LHS method is adopted, which generates the random sample (position of beetle) uniformly within a specified range depending upon the actual range of end milling process parameters. Algorithm 2 explains the modified LHS procedure.

Algorithm 2:

Initial sample generation using modified LHS

figure b

Furthermore, the newly generated beetle decides the next move based on the smell concentration by creating the next promising position in the neighborhood of previous position using the exploration and exploitation approach. The directional move is determined by Eq. (20),

$$ \overrightarrow{b}=\frac{rnd\left(k,1\right)}{\left\Vert rnd\left(k,1\right)\right\Vert } $$
(20)

where rnd is a random function, and k is the input dimension of the beetle position. The exploration is performed on right (xr) or left (xl) using Eq. (21) or Eq. (22),

$$ {x}_r={x}^t+{d}^t\times \overrightarrow{b} $$
(21)
$$ {x}_l={x}^t-{d}^t\times \overrightarrow{\mathrm{b}} $$
(22)

where dt is the sensing length of antennae at time t, which implies the exploitation skill. Value of dt must enfold the solution space. This guides the algorithm to escape from being stuck at local optima and improves the convergence speed. The next position of the beetle is decided using Eq. (23),

$$ {x}^t={x}^{t-1}+{\delta}^t\overrightarrow{b}\mathit{\operatorname{sign}}\left(f\left({x}_r\right)-f\left({x}_l\right)\right) $$
(23)

where δt is the step size of exploration, which follows a decreasing function of t, and sign represents a sign function. Sensing length dt and step size δt are updated using Eq. (24) and Eq. (25),

$$ {d}^t=0.95{d}^{t-1}+0.01 $$
(24)
$$ {\delta}^t=0.95{\delta}^{t-1} $$
(25)

The proposed data-driven BRNN-BAS framework is depicted in Fig. 9. It shows that the BRNN model is used as a surrogate model to the BAS optimization algorithm, which can evaluate candidate solutions effectively. The flowchart of the proposed BRNN-BAS algorithm is portrayed in Fig. 10. The sequential algorithmic flows are portrayed with dotted and continuous lines.

Fig. 9
figure 9

BRNN-BAS framework for end milling process

Full size image
Fig. 10
figure 10

BRNN-BAS algorithmic flowchart

Full size image

4.3.3 Performance metric

Root mean square error (RMSE) is used as the performance metric for the BRNN. The RMSE is an improved metric, which accurately measures regression errors. If the model-produced output response is y and the target response is t and i is the index of experiment for machining processes, the RMSE is calculated using,

$$ RMSE=\frac{1}{N}\sqrt{\sum_i{\left({y}_i-{t}_i\right)}^2}\kern0.5em $$
(26)

4.4 Regression model-assisted BAS technique

Multiple regression model is obtained for all the five responses at 95% confidence level. Regression equations are depicted in Eqs. (27)–(31). These equations are used as the fitness functions for the regression-BAS algorithm. Therefore, in Fig. 10, the dotted part is replaced using these equations for the regression-assisted optimization technique. Table 4 shows the p values and R2 values for the regression analysis, which states that the MRR is mostly dependent on the DOC and Ra is mostly dependent on the TD.

Table 4 Regression models with p values and R2 values
Full size table
$$ \mathrm{MRR}=-18.6+3.60\ast \mathrm{TD}-0.00464\ast \mathrm{SS}+0.12\ast \mathrm{FR}+12.73\ast \mathrm{DOC} $$
(27)
$$ \mathrm{Ra}=-0.885+0.1090\ast \mathrm{TD}+0.000290\ast \mathrm{SS}-0.1036\ast \mathrm{FR}+0.0937\ast \mathrm{DOC} $$
(28)
$$ \mathrm{Fx}=0.79+0.0434\ast \mathrm{TD}+0.000047\ast \mathrm{SS}+0.004\ast \mathrm{FR}-0.329\ast \mathrm{DOC} $$
(29)
$$ \mathrm{Fy}=9.88-0.239\ast \mathrm{TD}-0.00190\ast \mathrm{SS}-0.608\ast \mathrm{FR}-0.132\ast \mathrm{DOC} $$
(30)
$$ \mathrm{Fz}=17.52+1.200\ast \mathrm{TD}-0.00572\ast \mathrm{SS}-1.254\ast \mathrm{FR}-2.11\ast \mathrm{DOC} $$
(31)

5 Results and discussions

5.1 GRA analysis

In this study, GRA is applied on the data obtained from Table 3. Analysis is presented in Table 5. Using GRG, the multi-response problem becomes a single-response problem of maximization type. Taguchi’s analysis is applied on the GRG scores. The main effect plot for the parameters is presented in Fig. 11 a and mean response table for the GRG is shown in Table 6 (a). The near-optimal combination of the process parameters is TD3-SS3-FR3-DOC4, i.e., the values of TD = 8 mm, SS = 2000 rpm, FR = 4 mm/s, and DOC = 2 mm respectively. The delta values are portrayed (the difference of highest and lowest GRG mean) as 0.1680 for TD, 0.0654 for SS, 0.1705 for FR, and 0.1667 for DOC. This further implies that FR is the most influential parameter followed by TD, DOC, and SS for the end milling process. ANOVA results are displayed in Table 7.

Table 5 Grey relational analysis on Taguchi’s L16 OA design
Full size table
Fig. 11
figure 11

Main effect plots for a GRG and b VIKOR index

Full size image
Table 6 Response table for (a) GRG means and (b) VIKOR index means
Full size table
Table 7 ANOVA results on GRG
Full size table

ANOVA depicts the significance of the input variables on the GRG scores. It is observed that the TD and FR are most significant to the GRG, followed by the DOC and SS. It could be concluded that the SS is least significant to the process. Figure 12 shows the response surfaces for each two of the process variables vs. the obtained GRG values. For the TD-SS combination, a better GRG could be obtained when the TD varies between 8 and 9 mm and SS lies in the range of 1500–1600 rpm. For the TD-FR combination, higher GRG values are attained within certain range of TD values (8–9 mm) and the FR values (4–4.5 mm/s). When the DOC varies in the range of 1.75–2 mm and the TD remains in the same range, the GRG score improves.

Fig. 12
figure 12

Response surfaces of GRG vs. combinations of each two parameters

Full size image

The SS does not show much influence for the highest GRG, whereas the FR-DOC significantly affects the GRG scores within certain range of the FR (4–4–5 mm/s) and the DOC (1.5–2 mm). From the above study, the FR is shown to be most effective, and the SS is shown to be least effective for the end milling process.

5.2 VIKOR analysis

Five steps of the VIKOR method are applied on the multi-response data of Table 3 and the VIKOR index scores are obtained and presented in Table 8. Taguchi’s analysis is applied to the VIKOR index values. Equal weights for each of the responses are considered as, w1 = w2 = w3 = w4 = w5 = 0.2. The main effect plot for the parameters is presented in Fig. 11b. The mean response table for the VIKOR indices is shown in Table 6 (b).

Table 8 VIKOR method based on Taguchi’s L16 OA design
Full size table

The near-optimal setting for process parameters is TD2-SS4-FR2-DOC1, i.e., TD = 7 mm, SS = 2250 rpm, FR = 3 mm/s, and DOC = 0.5 mm respectively. The delta values in Table 6 (b) are 0.2678 for TD, 0.2010 for SS, 0.3623 for FR, and 0.3128 for DOC. This further implies that the FR is the most influential parameter followed by the DOC, TD, and SS for the end milling process. ANOVA is applied on the VIKOR indices and portrayed in Table 9. It is observed that the FR has the highest contribution to the VIKOR indices followed by the DOC, TD, and SS. It could be concluded that the SS is least significant to the process. Figure 13 shows the response surfaces based on the combination of two process variables and the corresponding VIKOR index values. The VIKOR index behavior is opposite of the GRG. For the TD-SS combination, a lower VIKOR index could be obtained when the TD varies between 8 and 9 mm and SS lies in the range of 1500–1600 rpm. For the TD-FR combination, lowest VIKOR indices are attained within specific ranges of TD (8–9 mm) and FR (4.4.5 mm/s). When the DOC varies in the range of 1.75–2 mm and the TD remains in the same range, the VIKOR indices are better. The SS does not show much influence on the VIKOR indices, which resembles with GRA study, whereas the FR-DOC is less significant to the VIKOR indices unlike the GRG scores. Here also, the FR is shown to be most effective and the SS is shown to be least effective to the process.

Table 9 ANOVA results on VIKOR index
Full size table
Fig. 13
figure 13

Response surfaces of VIKOR index values vs. combinations of each two parameters

Full size image

5.3 BRNN-BAS and regression-BAS comparison

The proposed BRNN-BAS and regression-BAS algorithms are applied on the data in Table 3 for the purpose of surrogate training and testing. The algorithms are coded in the MATLAB. The BRNN parameters are considered as the learning rate = 0.1, error goal = 1e−5, and number of epochs = 100. The data is divided in 70:30 for training and testing. The performance curve and regression plots of the BRNN are depicted in Fig. 14.

Fig. 14
figure 14

BRNN training and regression plots for the end milling data

Full size image

The training results and regression curves depict very high R values with low MSE, which proves that the BRNN model can approximate the function accurately. The cross validation is not considered since the high R values are obtained and accurate responses are predicted. The parameters of the BAS algorithm are set as d0 = 0.001, δ0 = 0.8, number of iterations = 500. The convergence plots of the data-driven BRNN-BAS algorithm and the regression-BAS algorithm are shown in Fig. 15 and Fig. 16 respectively. Both the algorithms are coded considering the boundary conditions of the process parameters. The parameters are evaluated based on the following range, 6 ≤ TD ≤ 10, 1500 ≤ SS ≤ 2250, 2 ≤ FR ≤ 5, and 0.5 ≤ DOC ≤ 2. Out of these four parameters, TD, SS, and FR are treated as the integer variables and DOC is treated as the real valued variable. The unfit solutions, which fall outside of the boundary conditions, are discarded from further evaluation. The multi-response end milling problem considered in this study is scaled down using the weighted sum method. The weights are derived from the VIKOR method. It could be observed that the proposed BRNN-BAS achieves the near-optimal solutions promptly. The best solution is fbst (objective value = − 6.2868) with parameters (TD = 10, SS = 1646, FR = 4.66, DOC = 1.54)t = 269. This solution shows improved responses (MRR = 46.027, Ra = 0.1652, Fx = 0.7823, Fy = 2.826, Fz = 10.8195).

Fig. 15
figure 15

Convergence curve obtained by BRNN-BAS

Full size image
Fig. 16
figure 16

Convergence curve obtained by regression-BAS

Full size image

5.4 Confirmatory test

Confirmatory tests are conducted to find the corrections for the GRG, VIKOR index, and objective scores obtained by both the BAS variants. The results are portrayed in Table 10.

Table 10 Confirmatory test results for GRA, VIKOR, regression-BAS, and BRNN-BSA methods
Full size table

The GRG value becomes 0.7791, which is close to the estimated value 0.8826. The VIKOR index is 0.8687, which is better than the predicted value of 1.2509. For the regression-BAS, the experimental objective score (− 3.4346) is inferior to the predicted score (− 4.7018). For the BRNN-BAS technique, the experimental run obtains an objective score of − 7.172, which is better than the predicted value − 6.2868. The responses obtained using the BRNN-BAS algorithm is (MRR = 42.042, Ra = 0.4651, Fx = 0.7973, Fy = 2.547, Fz = 2.3727). BRNN-BAS is shown to achieve 40.98% improved MRR and 10.56% improved Ra scores than the other methods. For the cutting forces, the obtained results are under satisfactory range. From the above study, it could be stated that the BRNN-BAS produces higher material removal rate, lower surface roughness, and the near-optimal set of cutting forces for the end milling on the Proxxon CNC machine.

6 Conclusions

In this article, a small desktop CNC milling machine is considered for the machining of the aluminum alloy used for sheet metal work. To obtain high-quality machining, a multi-response decision model is developed using the TD, SS, FR, and DOC as the process parameters. The MRR, Ra, and cutting forces (Fx, Fy, and Fz) are considered as the performance characteristics. Taguchi’s OAD coupled with GRA and VIKOR method, regression-assisted BAS algorithm, and BRNN-assisted BAS algorithm are considered for the process optimization task. Following conclusions are drawn from this study,

  • AA3105 shows substantially good machining properties under the end milling on the Proxxon FF 500 BL CNC with proper settings of the parameters. The ideal settings for the parameters are found to be (TD = 10 mm, SS = 1646 rpm, FR = 4.66 mm/s, and DOC = 1.54 mm).

  • The ANOVA and response surface plots depict that the FR is the most influential and SS is the least influential process parameters.

  • The improved GRG, VIKOR index, and objective values for both the BAS algorithms could be extremely helpful in minimizing the cost of manufacturing, which could improve the machining quality.

  • Probabilistic BRNN shows good performance curve and regression fit, which is comparable with the other predictive functions, such as the MLP, RBF, SVM, and Gaussian kernels. The BAS algorithm portrays its capability to obtain the near-optimal solutions promptly while coupled with the trained surrogate.

  • The BRNN-BAS technique shows its superior ability to optimize the machining process while compared with the heavily exploited GRA and VIKOR methods and regression-driven BAS algorithm. The BRNN-BAS is shown to outperform all the techniques by obtaining 40.98% improved MRR and 10.56% improved Ra values.

This work is aimed to extend with hard-to-machine materials, such as the steel alloys and composites as workpieces in future. More process variables are to be considered, such as radial depth of cut, axial depth of cut, tool wear, chatter vibration, power consumption by machines, and temperature. The many-objective optimization model for BAS would be developed for the said process, and comparisons would be performed with existing algorithms such as NSGA III or MOEA/D in future.