A novel MOALO-MODA ensemble approach for multi-objective optimization of machining parameters for metal matrix composites

Machining of metal matrix composites (MMCs) has now become an essential task in shaping them to their final usable forms for various industrial applications. These machining operations may range from drilling of simple through holes to cutting and shaping of the composites into complex shapes. Determination of the optimal process parameters during their machining presents a combinatorial optimization problem, which may turn out to be more complex due to involvement of various material dependent parameters of the MMCs, like particle size, reinforcement percentage etc. The importance of optimization increases manifold due to conflicting settings of different process parameters for attaining the desired quality characteristics. In this paper, two nature-inspired optimization algorithms, i.e. multi-objective antlion optimization (MOALO) and multi-objective dragonfly algorithm (MODA) are applied for optimization of various process parameters during machining of MMCs. To mitigate uncertainties in global optimality of the predicted Pareto optimal fronts arising due to stochastic nature of the metaheuristics, an ensemble approach integrating MOALO and MODA techniques is proposed here. It is demonstrated with the help of two case studies that MOALO-MODA ensemble is far superior to its individual counterparts and thus, can be employed for development of robust and reliable Pareto optimal fronts.

there are many parameters, like drill speed, feed rate, temperature, tool type, particle fraction etc. affecting the drilled surface quality and tool wear rate. Thus, the most influential parameters affecting efficiency of the machining operation and quality of the final products must be carefully determined. Basavarajappa et al. (2008) identified that during drilling of MMCs, the major factor influencing drilling force had been feed rate rather than cutting speed. It is also reported in the literature that surface roughness of the drilled holes on MMCs would increase with increase in feed rate and decrease with increase in cutting speed (Rubi et al. 2022;Davim and Monteiro Baptista 2001). But, many other researchers are against these observations (Tosun and Muratoglu 2004a;Ahamed et al. 2010). However, most of the researchers have agreed that feed rate would be the dominant parameter over cutting speed with respect to achieved surface roughness. Cutting tool material also plays an important role in influencing the drilled surface quality. For example, carbide coated tools would provide better surface quality than uncoated cutting tools (Basavarajappa et al. 2008) and HSS would not show any promising effect on drilled surface quality of MMCs. The HSS tools would often undergo excessive tool wear, resulting in rough surface during MMC drilling (Tosun and Muratoglu 2004b).
Presence of reinforcement has significant influence on the machining characteristics affecting machinability property of composites. Reinforcement volume fraction hinders chip formation due to higher plastic shear. During machining of MMCs, the hard ceramic particles when come in contact with the tool, start forming built-up edge (BUE). These small particles form a sandwich layer between the tool and the material, and the BUE acts as cutting tool causing generation of rough surface with higher tool wear. Machinability of MMCs can be effectively controlled while varying the reinforcement shape, size, properties of reinforcement, manufacturing method etc. (Ozben et al. 2008). As the matrix-embedded reinforcement particles are very hard and abrasive, machining of MMCs has become a challenging task with respect to excessive tool wear along with poor surface quality of the machined components. MMCs reinforced with silicon carbide particles are widely employed in industrial practice because of their easier and cheaper ex situ preparation process (Pramanik et al. 2007). Several investigations have been focussed on studying the machinability properties of MMCs reinforced with silicon carbide particles, mainly with respect to tool wear (Bhushan et al. 2010), surface integrity (Wang et al. 2022a) and chip formation (Guo et al. 2022).
Optimization of the machining processes is a crucial step in ensuring that MMCs can be efficiently and sustainably machined to their desired shapes (Nicholls et al. 2017). Determination of the optimal values of different machining parameters, like feed rate, depth of cut and speed helps in attaining better machinability without changing the properties of the considered materials. As surface quality of the machined composites plays a pivotal role in product performance, the concerned manufacturers should strictly attempt to adhere to its target value. Typically, surface roughness is considered as one of the main factors influencing surface integrity of the machined components (Outeiro 2014) According to Manna and Bhattacharayya (Manna and Bhattacharayya 2003), higher speed, lower depth of cut and lower feed rate would provide better surface finish during machining of Al/SiC composite. Similar optimal process parameter combinations have also been reported by the other researchers.
A significant number of studies has relied on the applications of various multi-criteria decision making (MCDM) tools for determination of the optimal intermixes of different process parameters during machining of MMCs. Sidhu and Yazdani (Sidhu and Yazdani 2018) employed desirability function and lexicographic goal programming to identify the optimal machining condition for having minimum residual stress and tool wear rate (TWR), and maximum material removal rate (MRR) during electrical discharge machining (EDM) of SiC/A359 MMC. Tamiloli et al. (2019) integrated fuzzy logic with technique for order of preference by similarity to ideal solution (TOPSIS) for optimization of end milling parameters while machining Al-SiC-fly ash MMC. It was observed that feed rate, speed and depth of cut would significantly affect surface roughness and cutting force during the said machining operation. Patel and Pradhan (2019) investigated the effects of various EDM process parameters on MRR, TWR, surface roughness and radial overcut during machining of Al7075/SiC/WS 2 composite. While machining AlSiC MMC, Satpathy et al. (2017) proposed the application of TOPSIS to find out the optimal levels of various EDM process parameters. Kumar et al. (2019) adopted additive ratio assessment (ARAS) method for optimality evaluation of EDM process parameters during machining of AA7050 composite. Kotlapati and Hiremath (2022) optimized EDM operation of Al-SiC MMC using grey relational analysis (GRA) considering hole depth, MRR, TWR, surface roughness and circularity error as the main responses.
Despite their significant advantages, like extremely low computational effort and ease of implementation, MCDM techniques have a major drawback that they rarely can find out the global optimal solution. They work best as a suitable alternative selection approach. On the other hand, several researchers have solved the machining process optimization problems treating them as single objective optimization (SOO) problems where only one response can be optimized. Senthilkumar et al. (2021) adopted response surface methodology to develop the corresponding metamodels for kerf width and surface roughness. Kumar et al. (2016a) applied genetic algorithm (GA) for optimizing the wear behaviour of Al2618 alloy reinforced with Si 3 N 4 , AlN and ZrB 2 particles. Different tribological properties, like wear rate, wear resistance, specific wear rate and coefficient of friction were considered to examine wear behaviour of the composite materials. Radhika et al. (2015) examined the wear behaviour of aluminium hybrid MMC and optimized the considered process parameters using GA. It was noticed that minimum wear rate for alumina composite could be achieved at 5 wt% reinforcement. Rao et al. (2021) applied GA and Taguchi methodology to optimize the drilling parameters for Al-SiC MMC. Bongale et al. (2018) attempted to minimize wear on Al/SiCnp/E-glass fibre MMC while adopting Taguchi methodology and GA technique. Adithiyaa et al. (2020) employed grey wolf optimizer and K-nearest neighbour algorithm for optimizing the process parameters of stir squeeze casting of AA2219 reinforced MMC. Kumar et al. (2022) proposed the application of TOPSIS to convert multiple objectives, like MRR, surface roughness and spark gap into a single composite objective function, which was later optimized using GA.
Although SOO is quite relevant for some applications where the main focus is on a particular response, in majority of the industrial applications, simultaneous achievement of multiple responses is often required. Thus, multi-objective optimization (MOO) needs special emphasis. However, many researchers have adopted formulations, like weighted sum approach (Bazgan et al. 2022;Chagas and Wagner 2022) to transform an MOO problem into a SOO problem. This is an efficient way to reduce the complexity and computational intensiveness of an optimization problem. However, among other challenges, specifying suitable weights to the responses a priori during the optimization phase is a major hurdle. The MOO generally predicts a set of optimal solutions known as non-dominated solutions which together form the Pareto optimal front. Such MOO involves two or more solutions that need to simultaneously maximized or minimized. Selvakumar et al. (2016) applied non-dominated sorting genetic algorithm-II (NSGA-II) to optimize the corresponding machining parameters that would minimize flank wear loss and maximize MRR of Al-4%Cu-7.5%SiC MMC. Singh et al. (2017) applied NSGA-II to optimize MRR and TWR during machining of MMCs. Kumar et al. (Kumar et al. 2016b) also employed NSGA-II to maximize MRR and minimize surface roughness while machining Al6061 silver-coated copper MMCs. Daniel et al. (2020) maximized MRR and minimized tool temperature during drilling operations of MMCs. Ikhlas et al. (2018) compared the predictive performance of artificial neural network and response surface methodology for metamodeling during machining operations of MMCs. The NSGA-II technique was also deployed to optimize surface roughness and cutting speed during machining of AISI 4140 hardened steel with mixed ceramic reinforcement particles. It was proved that NSGA-II would be more efficient than desirability function approach. Wu et al. (2022) endeavoured to combine deep learning with TOPSIS and GA to generate the corresponding Pareto fronts while solving machining optimization problems.
It can be revealed from the literature review of MOO applications in machining of MMCs that NSGA-II is the most commonly used algorithm. In this paper, two recently developed nature-inspired optimization algorithms, i.e. multi-objective antlion optimization (MOALO) and multi-objective dragonfly algorithm (MODA) are employed for multi-objective optimization of process parameters during machining of two different MMCs. The MOALO and MODA are then robustified while using an ensemble approach. A combined MOALO-MODA ensemble is finally proposed for the Pareto optimal front generation. The rest of the paper is organised as follows: Sect. 2 outlines the proposed ensemble approach of Pareto optimization. Discussions regarding MOALO, MODA and non-dominated sorting (NDS) are presented in Sects. 2.1, 2.2 and 2.3, respectively. Two case studies on machining of MMCs from the existent literature are solved using the proposed ensemble approach and are presented in Sect. 3. Section 3.1 demonstrates the effectiveness of this approach in predicting 2D Pareto optimal fronts for MOO during turning of MMCs. Similarly, in Sect. 3.2, 3D Pareto optimal fronts for MOO during drilling of MMCs are presented. Conclusions based on this comprehensive study are drawn in Sect. 4.

Methodology
In this paper, MOALO and MODA are applied for MOO of machining process parameters for MMCs. Taking inspiration from ensemble machine learning (ML) models, where a series of weak learners is considered to build a strong ML system, an ensemble approach to Pareto optimal front generation is proposed here. It is well known that to account for the stochastic nature of metaheuristic algorithms, multiple independent trials (n) of the algorithms must be carried out. Generally, the best Pareto optimal front with respect to spread and generational distance among n trials is reported. However, there may be other Pareto optimal fronts from those n trials that may have better non-dominated solutions at certain sections of the Pareto optimal front. To exploit this, an ensemble approach can be a suitable solution. A NDS algorithm based on NSGA-II is applied for secondary sorting in MOALO and MODA ensembles. Similarly, a tertiary NDS is also carried out in MOALO-MODA ensemble. It should be noted here that the primary NDS is carried out within each independent trial of MOALO and MODA.
The flowchart of the proposed methodology is illustrated in Fig. 1. The optimization process starts with specifying the objective function, design variables and their respective bounds. Then n independent optimization trials are carried out by both MOALO and MODA individually. The corresponding Pareto optimal fronts are predicted by MOALO and MODA in each of the n trials. The n Pareto optimal fronts predicted by MOALO are collated together and a NDS is carried out. The resulting Pareto optimal front is hereafter referred to as MOALO ensemble. Similarly, the Pareto optimal front obtained after NDS of n Pareto optimal fronts generated by MODA is referred to as MODA ensemble. The MOALO ensemble and MODA ensemble may be categorized as homogenous ensembles as the individual n Pareto optimal fronts that contribute to the final output are obtained using the same type of algorithm. Next, the ensemble Pareto optimal fronts of MOALO and MODA are collated together and a tertiary NDS is carried out to obtain MOALO-MODA ensemble Pareto optimal front. This ensemble is heterogeneous as it contains non-dominated solutions predicted by two different algorithms.

Multi-objective antlion optimization
Antlion optimizer (ALO) is a popular nature-inspired population-based metaheuristic algorithm, proposed by Mirjalili (2015) in 2015. The MOALO, the multi-objective version of ALO, was developed in 2017 by Mirjalili et al. (2017). It is based on the hunting tactics of an insect called the antlion. Antlions make traps for catching their prey by digging a cone-shaped pit in sand with their jaws. After making the pit, the antlion hides itself at the bottom of the cone and waits for the prey (primarily ants) to fall into the trap. Once the prey enters the pit, the antlion tries to consume it. However, the prey tries to evade the antlion and escape from the pit by crawling up. At the same time, the antlion starts tossing sand from the centre of the hole towards the edges so that the prey may slide back. Once the prey gets slided to the centre of the pit, it is unable to come out and antlion consumes it. To simulate the random walk of ants (prey), the following equation is used.
where cs calculates the cumulative sum, t shows the step of random walk and m is the maximum number of iterations.
To keep the random walk within the boundaries of the search space and letting the ant from overshooting the bounds, the following equation is considered for normalization.
where a i and b i are the minimum and maximum of random walks of ith variable, and c t i and d t i are the minimum and maximum of random walks of ith variable at tth iteration. The ALO simulates the entrapment of ants in antlion's pits by changing the random walks around antlion pits.
where c t and d t are the minimum and maximum of all variables at tth iteration, c t i and d t i are the minimum and maximum of all variables for ith ant, and A t j is the position of jth antlion at tth iteration. For simulating sliding of prey towards the antlion, the boundaries of the random walks are adaptively decreased.
where t and T are the current and maximum number of iterations, respectively, and w is the parameter responsible for adjusting accuracy of the exploitation. Based on the recommendations of Mirjalili et al. (2017), in this paper, w 2, 3, 4, 5 and 6, for (t/T ) > 0.1, 0.5, 0.75, 0.9 and 0.95 are, respectively, considered. Once the antlion catches the prey and consumes it, it reconstructs the pit.
where ant t i is the position of ith ant at tth iteration. Finally, elitism is introduced in ALO, wherein the fittest antlion achieved so far is recorded. The elite antlion influences the random walks of all the ants. Any ant performing random walk around an antlion which is selected by roulette wheel mechanism and the elite antlion can be mathematically modelled as below: where R t A and R t E are the random walks at tth iteration around the roulette wheel selected and elite antlions, respectively. An archive is used in MOALO to store the Pareto optimal solutions. To improve the Pareto optimal solution distribution and diversity in the archive, niching is adopted. The Pareto optimal solutions from the least populated neighbourhoods are chosen with the following probability: where c is a constant (greater than 1) and N i is the number of solutions near ith solution.
To accommodate new Pareto optimal solutions in the archive, in case if it is full, some solutions from the most populated regions are removed with the following probability: Further, for multi-objective optimization, Eq. (8) is modified as follows: Similarly, Eq. (9) is also modified for selection of antlions and elite. A roulette wheel mechanism and Eq. (9) are finally considered for selecting a Pareto optimal solution from the archive.

Multi-objective dragonfly algorithm
The MODA was introduced by Mirjalili (2016) in early 2015. This is also a nature-inspired population-based metaheuristic algorithm based on the behavioural patterns of dragonflies, mainly the hunting technique for food and swarming behaviour to save themselves from the enemy. Thus, it draws inspiration from the static and dynamic behaviours of the dragonflies. In static behaviour, they form small swarms in which they fly over a small area for hunting of other flying preys, like mosquitoes and butterflies. In dynamic behaviour, they form a very large swarm in which they migrate long distances in a particular direction. These static and dynamic swarming behaviours of the dragonflies are akin to exploration and exploitation, respectively. Movement of the dragonflies is governed by five main principles, i.e. separation, alignment, cohesion, attraction and distraction. Separation is the mechanism by which dragonflies avoid collision with the neighbouring dragonflies. It is mathematically modelled as: where S i is the separation of ith dragonfly, X i and X j are the current dragonfly and jth neighbouring dragonfly, respectively, and N is the number of dragonflies in the neighbourhood. Alignment is the mechanism by which dragonflies match their velocity to the neighbouring individuals, which can be modelled as below: ( 1 4 ) where A i is the alignment of ith dragonfly and V j is the velocity of jth neighbouring dragonfly. Cohesion is the mechanism by which dragonflies are attracted towards the centre of mass of the neighbourhood. It is mathematically modelled as follows: where C i is the cohesion of ith dragonfly. Attraction is the mechanism by which the dragonflies are attracted towards a food source. It can be mathematically expressed as below: where F i is the food source of ith dragonfly and X + is the location of the food source. Distraction is the mechanism by which the dragonflies move away from an enemy. It is mathematically modelled as: where E i is the enemy of ith dragonfly and Xis the location of the enemy. Based on these five behavioural patterns of the artificial dragonfly, step vector ( X) and position vector (X), the following relation is modelled to update its position.
where s a, c, f , e and w are the separation weight, alignment weight, cohesion weight, food factor, enemy factor and inertia weight, respectively, and t is the iteration counter. Next, the position vectors are updated using the following relation: In absence of any neighbourhood solution, dragonflies perform a Levy flight (Yang 2010) to update their positions.
where d is the dimension of the position vector. The Levy flights are computed as shown below: where r 1 and r 2 are random numbers between [0,1], β is a constant and σ is defined as follows: where Γ (x) (x -1)! An archive is used in MODA to store the Pareto optimal solutions. The position updating of search agents is identical to that of single objective DA, but the food sources are selected from the archive. To find out the least populated area of the Pareto optimal front, the search space is segmented. To achieve this, the best and the worst objective functions of the Pareto optimal solutions found so far are identified. A hyper-sphere is defined to cover all these solutions and the hyper-spheres are then divided into equal sub-hyperspheres in each iteration. The Pareto optimal solutions from the least populated neighbourhoods are chosen with the probability defined in Eq. (10). For selecting enemies from the archive, the worst (most populated) hyper-sphere is identified so that the artificial dragonflies are discouraged from searching around the non-promising crowded areas. Thus, some worst solutions from the most populated regions are removed with the probability defined in Eq. (11). To have an unbiased comparison between MOALO and MODA, the number of search agents for both of them is kept as 100, whereas, the iteration limit and primary archive size are kept as 100 and 500, respectively.

Non-dominated sorting
The NDS carried out in this paper during the secondary and tertiary sortings is the same as in NSGA-II (Deb et al. 2002). For secondary sorting, the Pareto optimal fronts obtained in n independent trials of MOALO (or MODA) are collated together and the number of solutions dominating each solution in the collated set is found out. Any solution B is said to be dominated by solution A if any the following three conditions is true (Sharma et al. 2017 The number of solutions (n d ) dominating each solution in the collated set is calculated. Any solution with n d 0 is assigned rank one. Here, n d 0 means that they are non-dominated solutions and are not dominated by any other solution in the collated set. Next, these solutions assigned rank as one are excluded from the collated set and the n d 0 in the reduced collated set is assigned rank two. This ranking procedure is iteratively carried out until all the solutions in the collated set are ranked. All the solutions in the collated set are then sorted according to their Pareto ranks assigned in the previous steps and the Pareto rank of Nth solution is found out (where N is the a priori set of maximum number of non-dominated solutions required for the ensemble Pareto optimal front). This rank is denoted as ndc. For all the solutions with Pareto rank ndc, the crowding distance in the objective function space is calculated. All the solutions are sorted according to their Pareto ranks, and then sorted with Pareto rank ndc in order of decreasing crowding distance. The first N solutions form the ensemble Pareto optimal set.

Metrics for Pareto optimal front comparison
The Pareto optimal fronts are assessed using four different metrics, i.e. generational distance (GD), convergence metric (CM), spread (SP) and inverse generational distance (IGD). The GD evaluates the closeness of non-dominated solutions to the known or true Pareto-optimal front (Veldhuizen et al. 1998). In this paper, to establish the superiority of MOALO-MODA ensemble, it is treated as the true Pareto optimal front. The GDs of MOALO ensemble and MODA ensemble are then compared with respect to MOALO-MODA ensemble. The value of GD is calculated using the following equation: where NP is the number of solutions in MOALO ensemble Pareto optimal front (or MODA ensemble Pareto optimal front) and d i is the Euclidean distance of each of these solutions to the nearest solution in the true Pareto optimal front. Similarly, CM too evaluates the closeness of non-dominated solutions to the true Pareto-optimal front (Zhou et al. 2006).
On the other hand, SP measures the distribution of nondominated solutions. It was originally propounded by Deb et al. (2011) and later extended by Zhou et al. (2006) as shown below: where M is the number of objective functions, e m is the number of M extreme solutions on the true Pareto optimal front, (e m ,S) is the Euclidean distance between the extreme solution on the true Pareto optimal front and extreme nondominated solution on MOALO ensemble Pareto optimal front (or MODA ensemble Pareto optimal front) based on mth objective function, i is the Euclidean distance of ith nondominated solution (S i ) to its nearest non-dominated solution (S j ) on MOALO ensemble Pareto optimal front (or MODA ensemble Pareto optimal front) and is the mean of i for all non-dominated solutions. The value of IGD (Wang et al. 2022b) is calculated based on the following expression: where NT is the number of solutions in the true Pareto optimal front and d i is the Euclidean distance of ith non-dominated solution on the true Pareto optimal front to its nearest solution in MOALO ensemble Pareto optimal front (or MODA ensemble Pareto optimal front). Further, for ease of comparison, the non-dominated solutions (S i ) predicted by the three ensembles are normalized (Ŝ i ) using Eq. (27) and the Pareto optimal front comparison metrics presented in Eqs. (23)-(26) are calculated again in their normalized forms.
where S T min and S T max are the lower and upper extremum nondominated solutions on the true Pareto optimal front.

Case study 1: turning of aluminium MMC
To demonstrate the efficacy of the proposed ensemble approach, a couple of case studies on machining of MMCs are considered here. Prakash et al. (2020) employed an L 27 orthogonal array to design a set of experiments during turning of rock dust reinforced aluminium MMC. The effects of reinforcement weight % (W ), particle size (P S ) (µm), turning speed (S) (rpm), feed rate (F) (mm/rev) and depth of cut (D C ) (mm) on MRR (g/min) and Ra (average surface roughness) (µm) were studied. During those experiments, each of the process parameters was varied at three different levels.
The experimental values of the considered process parameters and responses are presented in supplementary Table A1. Prakash et al. (2020) employed GRA technique to find out the optimal combination of the process parameters.
It is well documented in the literature that MCDM techniques are often unable to locate the global optimal solution, which leads to the deployment of metaheuristic algorithms for the said purpose. In this case study, the objective is to simultaneously maximize MRR and minimize Ra. However, due to conflicting nature of the parametric settings needed to individually maximize MRR and minimize Ra, instead of a unique solution, a set of non-dominated solutions called the Pareto optimal front is preferred. The proposed MOALO ensemble, MODA ensemble and MOALO-MODA ensemble are subsequently applied to determine the best Pareto optimal front. However, for deployment in the optimization scenario, the information (regarding dependency of the responses on the process parameters) derived from the experimental data in Table A1 must be expressed in the form of metamodels. For this purpose, polynomial regression is used to express MRR and Ra as functions of the considered process parameters, i.e. W , P S , S, F and D C . Initially, a full quadratic metamodel is developed and then step-wise elimination method is employed to remove the statistically insignificant terms from the metamodels. These insignificant terms in the model are identified based on their corresponding p-values (> 0.1) from analysis of variance (ANOVA) results. Thus, the following reduced quadratic metamodels are derived for MRR and Ra.
The corresponding ANOVA results for the reduced quadratic MRR and Ra metamodels are presented in Table  1. It can be noticed from this table that for MRR metamodel except for P S and W , the p-values for all other terms are less than 0.1, indicating them as statistically significant. It should be noted that these first-order terms (P S and W in MRR and W in Ra) with p-values greater than 0.1 are not eliminated to maintain hierarchical integrity of the metamodels. The responses predicted by the above-developed metamodels are compared with the experimental observations in Fig. 2. It is observed from Fig. 2 that all the data points lie on or just near the identity (diagonal) line indicating extremely high predictive power of the metamodels. To further examine these metamodels, their externally studentized residuals are plotted as normal probability plots with 95% confidence interval (CI) bands in Fig. 3. For MRR, all the residuals lie within the CI bands indicating a normal distribution of the residuals. On the other hand, for Ra, only one residual is noticed to be marginally outside the lower CI band. The goodness of fit of the metamodels is measured using R 2 , adjusted R 2 and predicted R 2 values, as reported in Table 2. There are negligible differences between R 2 and adjusted R 2 values for both the responses, which indicate that no statistically insignificant terms are included in the metamodels. Excellent values of predicted R 2 (which is the leave-one-out cross validation error metric) (99.15% and 96.69% for MRR and Ra, respectively) further confirm the accuracy of the developed  These metamodels are now deployed for optimization of the turning process parameters while machining the considered aluminium MMCs. All the subsequent optimization tasks are performed on an ASUS ROG Strix G531GT-AL030T series Windows system with Intel(R) Core™ i7 9th Gen, CPU @ 2.59 GHz, L2 Cache Size 1.5 MB and 8 GB RAM with 4 GB NVIDIA GeForce GTX1650. Five independent trials of MOALO and MODA are carried out, and the resulting Pareto optimal fronts are exhibited in Fig. 4. A lot of variability is noticed in the five different MOALO Pareto optimal fronts (Fig. 4a). It is unveiled from the box plots in Fig. 4a that the spread of the Pareto optimal fronts is inconsistent across various trials. The average number of non-dominated solutions in MOALO Pareto optimal fronts is 436 with a standard deviation of 94. A comparatively low (average) number of solutions as compared to the maximum archive size of 500 (for primary NDS) and relatively high standard deviation among various trials indicate that the results of MOALO are too stochastic. The best Pareto optimal front with respect to spread is obtained in trial 5 of MOALO. The average computational time for MOALO is 107 s with a standard deviation of 49 s. The curvature noticed at the lower side of some of the Pareto optimal fronts suggests that special care should be taken in choosing the optimal solutions from the lower Ra region. The Pareto optimal fronts generated by five different trials of MODA (Fig. 4b) are very similar to each other. Small variation is noticed only in the lower end (lower Ra region) of the Pareto optimal fronts of MODA. The average and standard deviation of the number of non-dominated solutions in MODA Pareto optimal fronts are 492 and 8, respectively. Each MODA run takes on an average 107 s (± 56 s) to complete. Thus, so far it can be concluded that MODA is more robust as compared to MOALO though both of them have similar computational time requirements. However, the largest range (i.e. difference between two extremum values of the Pareto optimal front) among the 10 trials is observed in case of MOALO trial.
The ensemble MOALO Pareto optimal front (generated using the five MOALO trials) and the ensemble MODA Pareto optimal front (generated based on the five MODA trials) are presented in Fig. 5a, b, respectively. It is worthwhile to mention here that a sufficiently high maximum archive size (5000) is considered in this paper during the secondary and tertiary NDS stages. This ensures that only the dominated solutions in the collated set are eliminated. It can be revealed from Fig. 5 that 747 and 2217 non-dominated solutions are present in MOALO and MODA ensembles, respectively. However, MOALO seems to have better nondominated solutions at lower end of the Pareto optimal front as compared to MODA. Thus, it becomes clear that despite high robustness of MODA, the utility of MOALO in navigating certain complex regions of the objective function cannot be ignored. This also justifies the need for heterogeneous ensemble.
Thus, the two homogenous ensembles discussed above are further passed through NDS for tertiary sorting that would eventually lead to the development of MOALO-MODA ensemble Pareto optimal front, presented in Fig. 6. The MOALO-MODA ensemble is made up of 2387 nondominated solutions out of which only 174 (~7.3%) are contributed by MOALO. Thus, as also observed from this figure, MODA accounts for majority of the non-dominated solutions in the heterogenous ensemble, but MOALO contributes some of the most important solutions towards the lower end of the Pareto optimal front. This is further highlighted in Fig. 7 using parallel plots for the ensembles. The solutions represented by red and blue colours vividly contrast the superiority of MODA in the middle and upper portions of the Pareto optimal fronts while MOALO solutions dominate the lower portion of the Pareto optimal fronts. It is also important to point out here that the secondary and tertiary NDSs consume less than 5 s each. As discussed earlier, for calculation of various metrics, like GD, CM, SP and IGD, MOALO-MODA ensemble is considered as the true Pareto optimal front. The MOALO and MODA ensembles are quantitatively compared with the true Pareto optimal front using these metrics in Table 3. It is noticed that MODA has better performance as compared to MOALO.
It should also be stated here that any Pareto optimal front contains a set of optimal non-dominated solutions, and the superiority of one cannot be arbitrarily established over the others. Thus, each solution of the Pareto optimal front is a compromise solution between two or more conflicting objectives. On the other hand, industrial problems generally require a single solution adhering to certain pre-specified needs. Perhaps, this may be the reason of much popularity of MCDM techniques in such applications, where objective or criteria weights (i.e. importance quotients) are implied a  priori. Nevertheless, it is hypothesized here and subsequently shown that MCDM application to Pareto optimal front can lead to much better solutions as opposed to MCDM-based level averaging approach widely prevalent in the literature. Table 4 shows the comparison of the derived solutions with those of Prakash et al. (2020). As stated earlier, GRA technique was employed by Prakash et al. (2020), which would not require objective weights to be explicitly described, but would calculate their values based on the decision matrix.
Thus, in this paper, GRA is also employed to determine the optimal solutions from the ensemble Pareto optimal fronts. It can be noticed that the solutions by MOALO ensemble are 14% and 99% better for MRR and Ra, respectively, as compared to the process parameters reported in Prakash et al. (2020) and calculated in this paper using the developed metamodels. Similarly, MODA ensemble solutions are found to be 51% and 10% better for MRR and Ra, respectively. The prediction of MOALO-MODA ensemble is the same as that  (Prakash et al. 2020) c Negative values indicate that the current results are poorer than (Prakash et al. 2020) of MOALO ensemble. Since MOALO-MODA ensemble is entirely made up of the solutions from both MOALO and MODA ensembles, it is often likely that the same solution as reported by one of the homogenous ensembles may be achieved.
To further show the effectiveness of the proposed ensemble approach in generating tailor-made solutions based on different end requirements, some additional example solutions are also reported in Table 4. A priori weights of 25%, 50% and 75% are assigned to MRR in three different scenarios. The TOPSIS method is applied here to determine the optimal solutions from the Pareto optimal fronts. Description of TOPSIS is beyond the scope of this paper, and can be found in Wang et al. (2020), Singh et al. (2020, and Kalita et al. (2022). In all these scenarios, the current solutions almost outperform the results reported in Prakash et al. (2020). It is also interesting to observe that at lower MRR weights (i.e. 25% and 50%), the TOPSIS-based optimal MOALO-MODA Pareto solutions are those which are drawn from MOALO ensemble. A sensitivity analysis based on TOPSIS scores is subsequently carried out in 99 trials by varying w MRR (weight assigned to MRR) from 1 to 99% in steps of 1% (sum of weights allocated to MRR and Ra responses must always be one). It is noticed that in 59 trials, the best solution is attained from MOALO ensemble and in 43 trials, it is achieved from MODA ensemble. In three trials (i.e. w MRR 0.4, 0.41 and 0.42), the best solution is obtained in both the ensembles. This further indicates that despite MODA's overall dominance in generating better non-dominated solutions, MOALO can also navigate certain portions of the objective function space better than MODA. It further emphasizes the need for the proposed heterogenous Pareto ensemble approach in process optimization.

Case study 2: drilling of boron carbide reinforced MMCs
In this case study, drilling operation of boron carbide reinforced MMC is considered. Based on Taguchi's L 27 orthogonal array, Taşkesen and Kütükde (Taşkesen and Kütükde 2014) conducted 27 experiments to determine thrust force (F a ) (N), drilling torque (M c ) (N cm) and surface roughness (Ra) (µm) while treating B 4 C particle weight % (W ), feed rate (F) (mm/rev), spindle speed (S) (rev/min) and drill material (D) as the input drilling parameters. The experimental values of those process parameters and the measured responses are presented in supplementary Table A2. During the experimental runs, Taşkesen and Kütükde (2014) considered three different levels for each of the process parameters and finally applied GRA technique to search out the optimal parametric combination.
To carry out multi-objective optimization for simultaneous minimization of F a , M c and Ra using the proposed ensemble technique, the following metamodels for the three  responses are subsequently developed: Since the focus of this paper is on development of a novel ensemble approach for multi-objective optimization and for the sake of compactness of this paper, the steps and results regarding ANOVA and other visual analyses of the metamodels are not presented here. These metamodels are now employed as objective functions to predict 3D Pareto optimal fronts.
At first, multi-objective optimization is carried for simultaneous minimization of F a , M c and Ra using MOALO and MODA. Five independent trials of MOALO and MODA are carried out, and the corresponding results are reported in the form of violin plots in Fig. 8. The three-dimensional (due to three responses) Pareto optimal fronts are easier to compare for each dimension, algorithm-wise using these violin plots. From Fig. 8a, it can be noticed that the 5th trial of MOALO is significantly different than the rest for F a values. On the other hand, MODA has similar distribution for all the trials. It is interesting to observe from Fig. 8a that in the all the trials of MOALO, ample number of non-dominated solutions exist in the upper value region of F a , which is almost inexistent in MODA trials. Similar observation is also noticed from Fig. 8b for M c . The stark differences between the nondominated solutions of MOALO and MODA are perhaps most prominent in Fig. 8c. From the distribution of the Pareto optimal fronts across the three dimensions for each trial, it is also become evident that while a certain trial may have better values in a particular dimension, other trials may have better solutions in some other dimension. This justifies the current strategy of collating multiple trials of the algorithm and obtaining an ensemble Pareto optimal front that would amalgamate the critical information of multiple trials. It should be noted here that MODA (average 72 s, standard deviation 15 s) is marginally less computationally intensive as compared to MOALO (average 87 s, standard deviation 14.7 s). Further, in both MOALO and MODA, for all the trials, the maximum archive size of 500 non-dominated solutions is reached. Figure 9 exhibits the homogenous ensemble 3D Pareto optimal fronts derived using MOALO and MODA. There are 1386 and 1885 non-dominated solutions in MOALO ensemble and MODA ensemble, respectively, after secondary NDS. As stated above, since the maximum archive size has been reached for every trial of MOALO and MODA, it signifies that the combined solution pool (of 5 trials) before secondary NDS is 2500 each for MOALO and MODA. Thus, MOALO ensemble and MODA ensemble reduce the NDS to 55.44% and 75.4% of the combined solution pool. Figure 9 also shows that MOALO has better solution density in higher M c and lower Ra regions, while the cluster of solutions for MODA is much better than MOALO in the mid Ra region. These contrasts in the solutions of MOALO ensemble and MODA ensemble in 3D objective space clearly indicate that there is sufficient scope for improvement of results for both the algorithms. Thus, tertiary NDS is adopted by collating MOALO ensemble and MODA ensemble solutions to obtain MOALO-MODA ensemble, presented in Fig. 10. Among 2868 non-dominated solutions in MOALO-MODA ensemble, 1124 (i.e.~39%) solutions are from MOALO ensemble. This is a much better equitable distribution of nondominated solutions than that observed in the previous 2D Pareto optimal case study. The parallel plots in Fig. 11 depict the distribution of the Pareto optimal solutions in the three ensembles for different drilling parameters and responses. These plots further help in understanding the solution pattern of the two different algorithms used in this paper. The quantitative comparison of the two homogeneous ensembles with the true Pareto optimal front is presented in Table 5, which proves that MODA has moderately better performance than MOALO.
The optimal solutions reported in Taşkesen and Kütükde (2014) are compared with the current solutions in Table 6. Minor improvements are noticed in F a for all the three ensembles, whereas, for M c , the current results are poorer than (Taşkesen and Kütükde 2014). However, the current Ra solutions are at least 54% better than the literature (Taşkesen and Kütükde 2014). It should again be stated here that the solutions selected using GRA (or any other MCDM technique) from the ensemble Pareto optimal fronts are strictly dependent on the objective weights. Thus, the ensemble Pareto optimal fronts provide enough opportunities to the decision  (Taşkesen and Kütükde 2014) c Negative values indicate that the current results are poorer than (Taşkesen and Kütükde 2014) maker (or end user) to choose a plethora of alternative solutions. For sake of brevity, unlike in previous case study, no additional solutions for different weights are shown in Table  6.

Conclusions
Multiple independent trials of metaheuristic optimization are usually carried out to account for the stochastic nature of these algorithms. The best Pareto optimal front with respect to spread and/or generational distance is then selected as a suitable solution to the optimization problem. However, the selected Pareto optimal front may have some suboptimal solutions in certain regions of the objective function space. Thus, in this paper, a novel ensemble approach to multi-objective optimization is proposed while considering MOALO and MODA algorithms. It relies on multiple re-runs of the metaheuristic optimization and then combines the nondominated solutions from these independent runs to form an ensemble Pareto optimal front. This is performed by collating the individual Pareto optimal fronts generated by the same algorithm and carrying out a secondary NDS to obtain the ensemble MOALO or ensemble MODA. These homogenous Pareto optimal fronts are further enhanced by collating the previous two ensembles and carrying out a tertiary NDS to obtain the heterogeneous MOALO-MODA ensemble. Based on two comprehensive case studies for multi-objective optimization of machining of MMCs, it is revealed that MODA is more robust and faster than MOALO. The non-dominated solution density of MODA is also observed to be better. However, MOALO is able to navigate certain sections of the objective function to locate effective non-dominated solutions that MODA had totally ignored. Nevertheless, MOALO solutions are accounted for only 7.3% and 39% of the total solutions in 2D and 3D MOALO-MODA ensemble Pareto fronts obtained in case studies 1 and 2, respectively. The run times of secondary and tertiary sorting are only a fraction of the total computational cost of the metaheuristic simulations. Thus, it is proposed that at practically negligible additional computational cost, better Pareto optimal fronts can be obtained by the proposed ensemble method. The limitation of this paper is in sticking to only two metaheuristics and that too from the same class of nature-inspired algorithms. However, in principle, this methodology can be extended to include Pareto optimal solutions from several other optimization algorithms. Enhancement of the Pareto optimal solutions by the proposed ensemble methodology while encompassing algorithms from different classes, like metaphor-less, gradient-based etc. may perhaps lead to more interesting results.

Declarations
Conflict of interest There is no conflict of interest among the authors. All the authors have gone through the manuscript and agreed to submit this paper.