# Application of Improved Monarch Butterfly Optimization for Parameters’ Optimization

Authors: Shixin Huang; Kedao Zhang; Hongmei Li; Xiangjian Chen
Journal: Mathematical Problems in Engineering
 (2023)

---

## Abstract

The reasonable selection of cutting parameters in the machining process is of great significance to improve productivity, reduce production costs, and improve the quality of parts. However, due to the complexity of cutting parameter model optimization, most factories currently use experience or refer to relevant manuals to select the value of cutting parameters in production. In order to avoid and minimize abnormalities, they usually select more experienced and conservative values, and often do not select reasonable cutting parameters, which is not conducive to improving productivity, reducing production costs, and improving the quality of parts. Therefore, the research on cutting parameter optimization has important theoretical value and application value. In this paper, in order to find the optimal cutting parameters, the cutting model is solved by the improved monarch butterfly optimization (IMBO) algorithm, and the optimized cutting parameters are obtained. By establishing the mathematical model of cutting, the constraint conditions of actual machining are introduced into the model. In order to solve the model, some ideas of particle swarm optimization (PSO) and differential evolution (DE) are added to the traditional monarch butterfly optimization (MBO) algorithm. The MBO algorithm is improved to deal with multiobjective optimization problems. The IMBO algorithm is used to optimize the cutting model. The experiment shows that the optimized cutting parameters can significantly reduce production cost and maintain high production efficiency. Compared with NSGA-II algorithm and other swarm intelligence optimization algorithms, it shows that the IMBO algorithm has certain advantages in multiobjective optimization.


---

## Body

## 1. Introduction

In the process of cutting, choosing appropriate cutting parameters can effectively reduce the production cost and improve the production efficiency. At present, enterprises usually choose the parameters according to the processing experience or the reference manual. These parameters are often conservative and cannot reach the optimal value. Although they can ensure that the parameters of processed parts meet the standards, they are not helpful to improve the production efficiency and reduce the production cost.With the rapid development of computer technology, the solution of cutting parameter model is also expanding rapidly. Swarm intelligence optimization algorithm has become an important tool to study cutting parameter optimization. At present, genetic algorithm and particle swarm optimization algorithm are widely used in China. Qi and Cheng of Northeast University [1], respectively, studied the optimization of milling and turning parameters based on genetic algorithms, and Ling et al. of Huazhong University of Science and Technology [2] studied the optimization of integrated parameters of multiple cutting methods based on genetic algorithms; Zhang and Chen of Beijing University of Aeronautics and Astronautics [3] and Liu and Huang of Tongji University [4], respectively, studied the optimization of cutting parameters based on particle swarm optimization, taking milling and turning as examples. Qin and Zhi of Wuhan University of Engineering [5] combined genetic algorithm and particle swarm optimization algorithm, the two optimization algorithms of biomimetics, to study cutting parameter optimization. Gao [6] of Huazhong University of Science and Technology compared various optimization algorithms applied to cutting parameter optimization problems, indicating that the hybrid particle swarm optimization algorithm with local search can obtain better results than the basic particle swarm optimization algorithm, genetic algorithm, etc. These studies can be summarized as follows: first, according to the expected goal of actual production, the mathematical model of cutting parameter optimization is established. Secondly, an optimization tool is sought to solve the established model. At present, swarm intelligence optimization algorithm has become the mainstream tool to solve the cutting parameter model. New intelligent optimization algorithms are constantly emerging. Each algorithm has its own advantages and disadvantages. How effective it is when applied to the optimization of cutting parameters needs to be compared and verified. Thirdly, because the cutting parameter optimization problem contains various complex nonlinear constraints, it is necessary to find a constraint processing mechanism. In recent years, swarm intelligence is used to solve optimization problems more and more frequently.In [7], a multiobjective particle swarm optimization algorithm based on two-archive mechanism (MOPSO_TA) was proposed for performing well in terms of convergence and diversity simultaneously. In [8], the paper presents a novel multiobjective particle swarm optimization feature selection method. In [9], a novel dynamic multiobjective optimization evolutionary algorithm based on particle swarm prediction strategy and prediction adjustment strategy was proposed. In [10], the author proposed an approach called genetically-modified multiobjective particle swarm optimization (GMPSO) for scheduling application workflows on hybrid clouds in the context of high-performance computing in an attempt to optimize Makespan and Cost. In [11], a hybrid particle swarm optimization algorithm, which employs an artificial fish swarm algorithm (AFSA) and a backward learning strategy in particle swarm optimization, is first proposed for constructing an optimal least squares support vector regression (OPT-LSSVR) model. In [12], this paper proposed a multisurrogate assisted binary particle swarm optimization, denoted as MS-assisted DBPSO. In [13], an optimization model aiming at processing time and energy consumption was established, and the PSO algorithm was used to solve the optimization problem. The effectiveness of the model was verified by experiments, but the constraints considered in the model were not comprehensive enough. Solarte-Pardo et al. [14] has taken the processing time as the goal and used the variable heterotopic genetic algorithm to solve the problem. It only considers the two constraints of the surface roughness of the part and the effective power of the machine tool, which is quite different from the actual processing constraints and does not consider the processing cost.Based on the abovementioned research results and existing problems, this article improves on the following two aspects.(1)
The cutting model was solved by the improved monarch butterfly optimization (IMBO) algorithm, and the optimized cutting parameters are obtained. By establishing the mathematical model of cutting, the constraint conditions of actual machining are introduced into the model.(2)
The model was solved using the IMBO algorithm. Some ideas for differential evolution algorithm and particle swarm optimization algorithm are added to accelerate the optimization speed and avoid falling into local optimum. The MBO algorithm is improved to deal with multiobjective optimization problems. The IMBO algorithm is used to optimize the cutting model.(3)
Consider the various constraints that exist in actual processing, use the constraint violation method to deal with constraints, and avoid the problem of the value of the penalty factor in the penalty function method from affecting the optimization results.The other chapters of this article are arranged as follows: a mathematical model of milling processing is established in Chapter 2, taking the lowest production cost and maximum production efficiency as the optimization goals, and the constraint conditions of actual processing are considered. The introductions of the traditional MBO algorithm and the IMBO algorithm are in Chapter 3. In Chapter 4, the IMBO algorithm is used to solve the milling parameter optimization model, and the optimized processing parameters are obtained and compared with other algorithms, it is proved that this algorithm has certain advantages compared with other algorithms. Chapter 5 is a summary of the full text.


## 2. Optimization of Milling Parameters

In the process of cutting, we should pursue lower production cost and higher processing efficiency on the premise of ensuring that the processing parameters meet the standard [14]. Reasonable selection of cutting parameters can effectively reduce production costs and improve production efficiency. At the same time, the selection of correct cutting parameters is also the guarantee of the quality of processed parts. However, due to very high complexity of the optimization of the cutting parameter model, the factory often sets the cutting parameters according to the machining manual or past experience during actual processing [15]. However, regardless of whether it is based on the machining manual or past experience, these cutting parameters are often conservative, which is of little help in reducing production costs and improving production efficiency. Therefore, the research on cutting parameter optimization has a very high benefit [16].


### 2.1. Cutting Parameter Optimization Model

In order to optimize the cutting parameters, the mathematical model of cutting parameters optimization should be established, including the design of variables, constraints, and objective function. The cutting parameters include cutting speedv, feed f, and cutting depth ap. Compared with the feed rate and cutting speed, the cutting depth has less influence on the system [17]. And the depth of cut can be determined by the remaining amount of workpiece machining and the specific requirements of actual machining. Therefore, the depth of cut can be treated as a known amount for subsequent processing.In this paper, the minimum production cost and maximum production efficiency are taken as the cutting parameter optimization objectives. During cutting, the total processing timeTw for one milling is calculated as follows [18]:(1)Tw=Tm+Th+Tc+Tot,where Tc is the tool change time of each process and Tot is the auxiliary time. Normally, these two times are constant and have little to do with the optimization of cutting parameters.Tm is the actual working time of the tool. The actual cutting time of the tool can be obtained by the following formula:(2)Tm=πDL1000vfZ,where Th is the tool change time of a process. The tool change time for a process is as follows:(3)Th=πLTR1000Cv1/mhv1/mhfyh/mh−1aePh/mhZuh/mhapkh/mhD1−mh/mh.In the above formulas,D is the diameter of the tool; L is the cutting length; v is the cutting speed; f is the feed per tooth of the milling cutter during milling; Z is the number of teeth of the milling cutter; TR is the tool change time; ae is the width of milling; ap is the depth of cut; and Cv, kh, yh, mh, uh, Ph, and qh are the wear resistance coefficients of the tool.At this time, the maximum production efficiency function is as follows:(4)minF1v,f=πDL1000vfZ+πLTR1000Cv1/mhv1/mhfyh/mh−1aePh/mhZuh/mhapkh/mhD1−mh/mh+Tc+Tot.When milling, the cost of a milling process is calculated as follows:(5)Cp=TwCtT0+C1+C0,where Ct is the cost of the tool, T0 is the wear resistance of the tool, C1 is the labor cost per unit time, and C0 is the management cost per unit time.Then, the minimum production cost function can be expressed as follows:(6)minF2v,f=TwCtT0+C1+C0.In production processing, due to the influence of processing materials, work piece quality, worker safety, spindle speed, and feed of the machine tool used, it is necessary to restrict the cutting speed, feed, cutting force, machine power, and surface roughness of parts.The cutting speed constraints are as follows:(7)g1v,f=πDwNmin60×1000−v≤0,g2v,f=v−πDwNmax60×1000≤0,where Dw is the diameter of the processed part, Nmin is the lowest spindle speed of the processing machine tool, and Nmax is the highest spindle speed of the processing machine tool.The feed rate constraints are as follows:(8)g3v,f=πDwvfmin1000Zv−f≤0,g4v,f=f−πDwvfmax1000Zv≤0,where vfmin is the minimum cutting feed speed of the machine tool and vmax is the maximum cutting feed speed of the machine tool.Cutting force constraint:(9)g5v,f=cFapxFfyFaeuFZNwFDwqFKFC−Ffmax≤0,where cF, xF, yF, uF, wF, qF, and KFC are the cutting force parameters and Ffmax is the maximum force of machine tool spindle.Machine tool power constraint:(10)g6v,f=kpvfap1000η−Pcmax≤0,where kp is the scale factor, Pcmax is the maximum effective cutting power, and η is the rotation coefficient of the machine tool.The constraint expression of the roughness requirement of the part surface is as follows:(11)g7v,f=1000f28rε−Ramax≤0,where rε is the arc radius of the tool tip and Ramax is the maximum surface roughness of the part.After the above constraint analysis, the optimization model of cutting parameters can be summarized as follows:(12)minFO1,O2,st.giv,f≤0i=1,2,…,7,where O1 is F1v,f and O2 is F2v,f.


### 2.2. Constraint Treatment

From the above established cutting parameter optimization model, it can be seen that there are a lot of nonlinear constraints in cutting parameter optimization. When using the IMBO algorithm to optimize cutting parameters, corresponding constraint processing is required. The common way to deal with constraints is to transform constrained optimization problems into equivalent unconstrained optimization problems. At present, constraint violation degree method [19] is widely used to deal with unconstrained optimization problems.The expression of constraint violation functiongx is as follows:(13)GX=∑i=1Pnmax0,giX+∑j=1qnmax0,hjX−δ,where X is the function variable, giX is the i-th inequality constraint, Pn is the number of inequality constraints, hjX is the j-th equality constraint, qn is the number of equality constraints, and δ is the shape coefficient, which is usually a real number with a very small value. Here, take δ=0.0001. The greater the value of the constraint violation function GX, the worse the feasibility of the variable X, and its magnitude is the degree of constraint violation.This paper will use the constraint violation degree method to deal with the constraint conditions.


## 2.1. Cutting Parameter Optimization Model

In order to optimize the cutting parameters, the mathematical model of cutting parameters optimization should be established, including the design of variables, constraints, and objective function. The cutting parameters include cutting speedv, feed f, and cutting depth ap. Compared with the feed rate and cutting speed, the cutting depth has less influence on the system [17]. And the depth of cut can be determined by the remaining amount of workpiece machining and the specific requirements of actual machining. Therefore, the depth of cut can be treated as a known amount for subsequent processing.In this paper, the minimum production cost and maximum production efficiency are taken as the cutting parameter optimization objectives. During cutting, the total processing timeTw for one milling is calculated as follows [18]:(1)Tw=Tm+Th+Tc+Tot,where Tc is the tool change time of each process and Tot is the auxiliary time. Normally, these two times are constant and have little to do with the optimization of cutting parameters.Tm is the actual working time of the tool. The actual cutting time of the tool can be obtained by the following formula:(2)Tm=πDL1000vfZ,where Th is the tool change time of a process. The tool change time for a process is as follows:(3)Th=πLTR1000Cv1/mhv1/mhfyh/mh−1aePh/mhZuh/mhapkh/mhD1−mh/mh.In the above formulas,D is the diameter of the tool; L is the cutting length; v is the cutting speed; f is the feed per tooth of the milling cutter during milling; Z is the number of teeth of the milling cutter; TR is the tool change time; ae is the width of milling; ap is the depth of cut; and Cv, kh, yh, mh, uh, Ph, and qh are the wear resistance coefficients of the tool.At this time, the maximum production efficiency function is as follows:(4)minF1v,f=πDL1000vfZ+πLTR1000Cv1/mhv1/mhfyh/mh−1aePh/mhZuh/mhapkh/mhD1−mh/mh+Tc+Tot.When milling, the cost of a milling process is calculated as follows:(5)Cp=TwCtT0+C1+C0,where Ct is the cost of the tool, T0 is the wear resistance of the tool, C1 is the labor cost per unit time, and C0 is the management cost per unit time.Then, the minimum production cost function can be expressed as follows:(6)minF2v,f=TwCtT0+C1+C0.In production processing, due to the influence of processing materials, work piece quality, worker safety, spindle speed, and feed of the machine tool used, it is necessary to restrict the cutting speed, feed, cutting force, machine power, and surface roughness of parts.The cutting speed constraints are as follows:(7)g1v,f=πDwNmin60×1000−v≤0,g2v,f=v−πDwNmax60×1000≤0,where Dw is the diameter of the processed part, Nmin is the lowest spindle speed of the processing machine tool, and Nmax is the highest spindle speed of the processing machine tool.The feed rate constraints are as follows:(8)g3v,f=πDwvfmin1000Zv−f≤0,g4v,f=f−πDwvfmax1000Zv≤0,where vfmin is the minimum cutting feed speed of the machine tool and vmax is the maximum cutting feed speed of the machine tool.Cutting force constraint:(9)g5v,f=cFapxFfyFaeuFZNwFDwqFKFC−Ffmax≤0,where cF, xF, yF, uF, wF, qF, and KFC are the cutting force parameters and Ffmax is the maximum force of machine tool spindle.Machine tool power constraint:(10)g6v,f=kpvfap1000η−Pcmax≤0,where kp is the scale factor, Pcmax is the maximum effective cutting power, and η is the rotation coefficient of the machine tool.The constraint expression of the roughness requirement of the part surface is as follows:(11)g7v,f=1000f28rε−Ramax≤0,where rε is the arc radius of the tool tip and Ramax is the maximum surface roughness of the part.After the above constraint analysis, the optimization model of cutting parameters can be summarized as follows:(12)minFO1,O2,st.giv,f≤0i=1,2,…,7,where O1 is F1v,f and O2 is F2v,f.


## 2.2. Constraint Treatment

From the above established cutting parameter optimization model, it can be seen that there are a lot of nonlinear constraints in cutting parameter optimization. When using the IMBO algorithm to optimize cutting parameters, corresponding constraint processing is required. The common way to deal with constraints is to transform constrained optimization problems into equivalent unconstrained optimization problems. At present, constraint violation degree method [19] is widely used to deal with unconstrained optimization problems.The expression of constraint violation functiongx is as follows:(13)GX=∑i=1Pnmax0,giX+∑j=1qnmax0,hjX−δ,where X is the function variable, giX is the i-th inequality constraint, Pn is the number of inequality constraints, hjX is the j-th equality constraint, qn is the number of equality constraints, and δ is the shape coefficient, which is usually a real number with a very small value. Here, take δ=0.0001. The greater the value of the constraint violation function GX, the worse the feasibility of the variable X, and its magnitude is the degree of constraint violation.This paper will use the constraint violation degree method to deal with the constraint conditions.


## 3. Optimization Algorithm

This chapter first introduces the traditional MBO algorithm. Then, some ideas of PSO and DE are added to the traditional MBO algorithm. Finally, the MBO algorithm is improved to deal with the multiobjective optimization problem.


### 3.1. MBO Algorithm

Monarch butterflies migrate from southern Canada and the northern United States to Mexico each summer for winter, traveling thousands of miles. The migration started in August and ended with the first frost. In spring they fly back from Mexico to the United States and Canada. Female butterflies lay eggs and produce offspring during migration. According to the migration behavior of monarch butterflies, Wang et al. [20] proposed the optimization algorithm of monarch butterflies.In the Monarch Butterfly Optimization algorithm, all butterflies are idealized and follow the following rules: (1) All butterflies can only be in southern Canada, northern United States (Land1) and Mexico (Land2), which means that the butterflies in these two areas are the whole butterfly population. (2) All butterfly offspring are generated by the butterflies in Land1 or Land2 after migration. (3) In order to ensure the size of butterfly population is fixed, it is stipulated that if the offspring butterflies are better in fitness than the parent generation, the parent generation shall be replaced by the offspring; if the offspring are less fit than the parent, the offspring are discarded directly. (4) The most suitable butterfly will be retained and directly enter the next generation to ensure that the population quality will not decline.The MBO algorithm simulates the migration process of monarch butterfly, and updates the candidate solution through the information exchange among butterflies. After simplification and idealization, all butterfly individuals in MBO are only in Land1 and Land2, and butterfly position updating is realized by migration and adjustment. In other words, the position of offspring butterflies is determined by migration operation, and other butterflies update their position by adjustment operation. Because the adjustment operation and migration operation can be executed simultaneously, the MBO can be processed in parallel, and the global and local search ability can be considered well.In MBO algorithm, after initializing all the parameters, the initial population should be generated and the fitness value of the initial population should be calculated; the butterflies are sorted according to the fitness value, and then the positions of the butterflies are updated according to the rules until the termination condition is met. In order to ensure the stability of population size and reduce the calculation times of fitness value as much as possible, the butterflies in Land1 and Land2 were moved and adjusted to update their positions.Migration operation:The migration operation is to update the position of the monarch butterflies between Land1 and Land2. If the total number of butterflies is set asNP and the butterfly mobility is p, then the number of butterflies in land1 is ceilp×NP, and the number of butterflies in land2 is NP−ceilp×NP. ceilx means to round x to an integer. For the convenience of expression, butterflies at Land1 and Land2 are Subset1 and Subset2. Then, subset1 = ceilp×NP and subset2 = NP−ceilp×NP.The migration of butterflies can be expressed as follows:(14)xi,kt+1=xr1,kt,r≤p,xr2,kt,r>p,where t represents the current number of iterations, r1 means to take a butterfly randomly from Subset1, r2 means to take a butterfly randomly from Subset 2, xi,kt+1 represents the k-th dimension of xi in the t+1 iteration, xr1,kt represents the k-th dimension of xr1 in t iteration, and xr2,kt represents the k-th dimension of xr2 in t iteration. r=rand×peri, where rand is a random number from 0 to 1. peri is the migration period, and in this case because there are 12 months in a year, it is usually 1.2.According to formula, the migration ratep determines the direction of butterfly population migration operation. If the migration rate is large, most of the newly generated butterflies appear in Subset1. Subset2 had more offsprings if the migration rate is low. According to the reality that butterflies have 5 months on Land1, set p to 5⁄12 [21].Adjustment action: the migration operation updates the parent with the offspring, while the adjustment operation updates the position of the butterfly individuals without generating new butterfly individuals.For individualj of a butterfly, a random number rand is generated. If rand≤p, butterfly j updates the position information according to the following formula:(15)xj,kt+1=xbest,kt,where xj,kt+1 represents the k-th dimension of xj in t+1 iteration, and xbest,kt is the k-th dimension of butterfly individuals in the optimal position in Land1 and Land2 in t iterationIfrand>p, butterfly j updates the position information according to the following formula:(16)xj,kt+1=xr3,kt,where xj,kt+1 represents the k-th dimension of xj in t+1 iteration. xr3,kt represents the k-th dimension of xr3 in t iteration. Butterfly r3 is randomly selected from Land2. At this point, if rand>BAR, you need to continue updating as follows:(17)xj,kt+1=xj,kt+1+α×dxk−0.5,where BAR is the adjustment rate, BAR=p. α is the weight factor, α=Smax/t2. Smax is the maximum step of a butterfly. dx is the step of butterfly j, dx=Levyxjt. The formula shows that the larger the α, the longer the search step and the stronger the global search ability. Conversely, the smaller the α, the smaller the search step and the better the local search ability.Basic algorithm flow: from the above analysis, it is known that in MBO algorithm, after initializing the position information of monarch butterfly, monarch butterfly individuals continuously update the position information through migration and adjustment operations until a certain condition is met. Here are the basic steps of the MBO algorithm:Step 1: Initialize. At this time, the number of iterationst=1, set the maximum number of iterations, and randomly generate NP monarch butterfly individuals. Set Smax, BAR, p, peri, and other parameters.Step 2: Calculate the fitness value. The fitness values are calculated according to the position information of the initialized monarch butterfly.Step 3: Divide the monarch butterfly population into two subpopulations.Step 4: Perform migration and tuning operations.Step 5: Merge the two subpopulations again.Step 6: Recalculate the fitness values of all monarch butterflies.Step 7:t=t+1, judge whether the stop condition is met, if so, the optimal result will be output; if not, skip to Step 3.


### 3.2. Fast Nondominant Sorting

According to [22], for any solution p, two properties are calculated: first, the number of solutions dominate the solution p, which is recorded as np; second, the set of solutions called Sp dominated by the solution p. The calculation steps are as follows:Step 1: Leti=0, the number of dominant solutions np=0. The set of dominant solutions Sp=∅.Step 2: For each individualq in the set P, if p≺q, individual q is added to the set Sp of solutions dominated by solution p; if q≺p, the number of solutions dominating p is np=np+1. Finally, we can get the dominating number np of each individual and the set Sp of dominating solutions. The individuals with np=0 are put into the front F1, so that the nondominated rank of p is Prank=1.Step 3: LetQ=∅, for each individual p in Fi, do the following: For individual q in Sp, nq=nq−1. If nq=0, then the nondominated level of Q is qrank=i+1, and the individual q is added to the set Q.Step 4: ifQ≠∅, i=i+1, Fi=Q, turn to Step 3; if Q=∅, stop.


### 3.3. Multiobjective MBO Algorithm

In this paper, the MBO algorithm is applied to milling parameter optimization, which involves multiple optimization objectives. The MBO introduced in the previous chapter is single-objective, so it needs to be transformed into a multiobjective. The modified solution steps are as follows:First, initialize. Randomly generate the parent populationP0 with t=0. The objective function values of each solution in the population P0 are calculated and assigned to the corresponding solution.Second, save the nondominated solutions found during the search to the external fileRep. There are two cases for preserving a nondominant solution at this point: if the external file Rep is empty, it directly determines whether the solution is a nondominant solution. If the solution is nondominant, it is added to the external file Rep. If the external file Rep is not empty, first determine if the solution is a nondominant solution, then determine if the solution is not dominated by the solution in the external file Rep. If it is satisfied, add the solution to the external file Rep, and then delete the other solutions that are dominated by the solution.Third, IMBO algorithm was used to generate the next generation population. Repeat the above process until the stop condition is met.For the second step above, the process of storing the nondominant solution is described in detail below. Figures1 and 2 show the storage of nondominant solutions when the external file Rep is empty and not empty, respectively.Figure 1
Storage when the external file is empty.Figure 2
Storage when the external file is not empty.Figure1 shows the storage situation when the external file Rep is empty. At this point, it is directly judged whether the solution is a nondominant solution. If it is not, it is directly discarded (as shown in Case 1). If the solution is nondominant, it is added directly to the external file Rep (as shown in Case 2). Figure 2 shows the case when the external file is not empty. In this case, the first step is to determine whether the solution is a nondominant one. If the solution is not a nondominant solution, it is discarded directly (as shown in Case 1). If the solution is a nondominant solution, it needs to be further determined whether it is dominated by an existing solution in the external file Rep. If the solution is dominated by the solution in an external file, it is discarded (as shown in Case 1). If the solution is not controlled by an existing solution in the external file, it is added to the external file Rep. At this time, delete the other solutions in the external file that are dominated by this solution. Case 2 in Figure 2 above indicates that there is no other solution governed by this solution in the external file, so it is added directly to the external file Rep. Case 3 indicates that x2 and x3 are dominated by the new solution after adding the new solution to the external file. In this case, x2 and x3 need to be removed from the external file.When the MBO algorithm is updated iteratively, the corresponding dimension values of randomly selected individuals are directly assigned to the new individuals, which limits the diversity of the population and is detrimental to its evolution. For this reason, some ideas of differential evolution and particle swarm optimization (PSO) are added to the MBO algorithm to increase the diversity of the population. In differential evolution algorithms, the values on the corresponding dimensions of two individuals are subtracted to produce a difference vector, which is assigned a weight and added to another vector to produce a variation vector, thus increasing the diversity of the population. PSO algorithms update iteratively according to the individual extremes and global extremes, which can speed up convergence. Based on the above analysis, add the following actions to the migration operation:(18)xj,kt+1=xpbest,kt+xneighbor1,kt−xneighbor2,kt×R.In the above formula,xj,kt+1 is the k-th dimension of xj at the t+1 iteration, xpbest,kt is the k-th dimension of the individual extremum at the t-th iteration and xneighbor1t and xneighbor2t are two different individuals randomly selected from Subset2 during the t-th iteration. R is the weight assigned to the difference vector and can be calculated from the following formula:(19)R=rand−0.5×2.As can be seen from the above formula,R is a random number between −1 and 1. The larger the R, the greater the influence of the difference vector on the variation vector.NSGA-II uses a fast nondominant sorting method, which reduces computational complexity. Adding some operations in NSGA-II to the MBO algorithm can achieve better performance.During the iteration, the populations are merged after the migration and adjustment operations. The combined population was placed in a nondominated order, and the crowding degree of each individual was calculated. Then, the parent individuals are selected according to the crowding degree, and the solution is updated by the crossover operator. The following formula is the calculation method of new individualsXA′ and XB′ after the arithmetic crossover operation of XA and XB:(20)XA′=aXB+1−aXA,XB′=aXA+1−aXB.In the formula above,a is a random number between −0.5 and 1.5. The above operation can better maintain the diversity of the population.The flow chart of IMBO algorithm is shown in Figure3. The specific implementation process is as follows: first, initialize the population, external file Rep, population size, and other parameters. Then, the fitness values of all individuals in the population are calculated, and the nondominated solutions are found and stored in Rep; record the best position of the individual pbest, and the global optimal solution gbest. Finally, the algorithm is executed iteratively. The population is divided into two subpopulations, namely subset1 and subset2, which are transferred and adjusted to generate new solutions. Two new subpopulations are merged and nondominated sorting is performed. The crowding degree is calculated at the same time. Then, the parent individuals are selected for arithmetic crossover operation to further update the solution. The fitness of each individual is calculated and Rep, pbest, and gbest are updated. When the termination condition is satisfied, the algorithm is terminated; otherwise, continue the iteration.Figure 3
Schematic flowchart of IMBO.


## 3.1. MBO Algorithm

Monarch butterflies migrate from southern Canada and the northern United States to Mexico each summer for winter, traveling thousands of miles. The migration started in August and ended with the first frost. In spring they fly back from Mexico to the United States and Canada. Female butterflies lay eggs and produce offspring during migration. According to the migration behavior of monarch butterflies, Wang et al. [20] proposed the optimization algorithm of monarch butterflies.In the Monarch Butterfly Optimization algorithm, all butterflies are idealized and follow the following rules: (1) All butterflies can only be in southern Canada, northern United States (Land1) and Mexico (Land2), which means that the butterflies in these two areas are the whole butterfly population. (2) All butterfly offspring are generated by the butterflies in Land1 or Land2 after migration. (3) In order to ensure the size of butterfly population is fixed, it is stipulated that if the offspring butterflies are better in fitness than the parent generation, the parent generation shall be replaced by the offspring; if the offspring are less fit than the parent, the offspring are discarded directly. (4) The most suitable butterfly will be retained and directly enter the next generation to ensure that the population quality will not decline.The MBO algorithm simulates the migration process of monarch butterfly, and updates the candidate solution through the information exchange among butterflies. After simplification and idealization, all butterfly individuals in MBO are only in Land1 and Land2, and butterfly position updating is realized by migration and adjustment. In other words, the position of offspring butterflies is determined by migration operation, and other butterflies update their position by adjustment operation. Because the adjustment operation and migration operation can be executed simultaneously, the MBO can be processed in parallel, and the global and local search ability can be considered well.In MBO algorithm, after initializing all the parameters, the initial population should be generated and the fitness value of the initial population should be calculated; the butterflies are sorted according to the fitness value, and then the positions of the butterflies are updated according to the rules until the termination condition is met. In order to ensure the stability of population size and reduce the calculation times of fitness value as much as possible, the butterflies in Land1 and Land2 were moved and adjusted to update their positions.Migration operation:The migration operation is to update the position of the monarch butterflies between Land1 and Land2. If the total number of butterflies is set asNP and the butterfly mobility is p, then the number of butterflies in land1 is ceilp×NP, and the number of butterflies in land2 is NP−ceilp×NP. ceilx means to round x to an integer. For the convenience of expression, butterflies at Land1 and Land2 are Subset1 and Subset2. Then, subset1 = ceilp×NP and subset2 = NP−ceilp×NP.The migration of butterflies can be expressed as follows:(14)xi,kt+1=xr1,kt,r≤p,xr2,kt,r>p,where t represents the current number of iterations, r1 means to take a butterfly randomly from Subset1, r2 means to take a butterfly randomly from Subset 2, xi,kt+1 represents the k-th dimension of xi in the t+1 iteration, xr1,kt represents the k-th dimension of xr1 in t iteration, and xr2,kt represents the k-th dimension of xr2 in t iteration. r=rand×peri, where rand is a random number from 0 to 1. peri is the migration period, and in this case because there are 12 months in a year, it is usually 1.2.According to formula, the migration ratep determines the direction of butterfly population migration operation. If the migration rate is large, most of the newly generated butterflies appear in Subset1. Subset2 had more offsprings if the migration rate is low. According to the reality that butterflies have 5 months on Land1, set p to 5⁄12 [21].Adjustment action: the migration operation updates the parent with the offspring, while the adjustment operation updates the position of the butterfly individuals without generating new butterfly individuals.For individualj of a butterfly, a random number rand is generated. If rand≤p, butterfly j updates the position information according to the following formula:(15)xj,kt+1=xbest,kt,where xj,kt+1 represents the k-th dimension of xj in t+1 iteration, and xbest,kt is the k-th dimension of butterfly individuals in the optimal position in Land1 and Land2 in t iterationIfrand>p, butterfly j updates the position information according to the following formula:(16)xj,kt+1=xr3,kt,where xj,kt+1 represents the k-th dimension of xj in t+1 iteration. xr3,kt represents the k-th dimension of xr3 in t iteration. Butterfly r3 is randomly selected from Land2. At this point, if rand>BAR, you need to continue updating as follows:(17)xj,kt+1=xj,kt+1+α×dxk−0.5,where BAR is the adjustment rate, BAR=p. α is the weight factor, α=Smax/t2. Smax is the maximum step of a butterfly. dx is the step of butterfly j, dx=Levyxjt. The formula shows that the larger the α, the longer the search step and the stronger the global search ability. Conversely, the smaller the α, the smaller the search step and the better the local search ability.Basic algorithm flow: from the above analysis, it is known that in MBO algorithm, after initializing the position information of monarch butterfly, monarch butterfly individuals continuously update the position information through migration and adjustment operations until a certain condition is met. Here are the basic steps of the MBO algorithm:Step 1: Initialize. At this time, the number of iterationst=1, set the maximum number of iterations, and randomly generate NP monarch butterfly individuals. Set Smax, BAR, p, peri, and other parameters.Step 2: Calculate the fitness value. The fitness values are calculated according to the position information of the initialized monarch butterfly.Step 3: Divide the monarch butterfly population into two subpopulations.Step 4: Perform migration and tuning operations.Step 5: Merge the two subpopulations again.Step 6: Recalculate the fitness values of all monarch butterflies.Step 7:t=t+1, judge whether the stop condition is met, if so, the optimal result will be output; if not, skip to Step 3.


## 3.2. Fast Nondominant Sorting

According to [22], for any solution p, two properties are calculated: first, the number of solutions dominate the solution p, which is recorded as np; second, the set of solutions called Sp dominated by the solution p. The calculation steps are as follows:Step 1: Leti=0, the number of dominant solutions np=0. The set of dominant solutions Sp=∅.Step 2: For each individualq in the set P, if p≺q, individual q is added to the set Sp of solutions dominated by solution p; if q≺p, the number of solutions dominating p is np=np+1. Finally, we can get the dominating number np of each individual and the set Sp of dominating solutions. The individuals with np=0 are put into the front F1, so that the nondominated rank of p is Prank=1.Step 3: LetQ=∅, for each individual p in Fi, do the following: For individual q in Sp, nq=nq−1. If nq=0, then the nondominated level of Q is qrank=i+1, and the individual q is added to the set Q.Step 4: ifQ≠∅, i=i+1, Fi=Q, turn to Step 3; if Q=∅, stop.


## 3.3. Multiobjective MBO Algorithm

In this paper, the MBO algorithm is applied to milling parameter optimization, which involves multiple optimization objectives. The MBO introduced in the previous chapter is single-objective, so it needs to be transformed into a multiobjective. The modified solution steps are as follows:First, initialize. Randomly generate the parent populationP0 with t=0. The objective function values of each solution in the population P0 are calculated and assigned to the corresponding solution.Second, save the nondominated solutions found during the search to the external fileRep. There are two cases for preserving a nondominant solution at this point: if the external file Rep is empty, it directly determines whether the solution is a nondominant solution. If the solution is nondominant, it is added to the external file Rep. If the external file Rep is not empty, first determine if the solution is a nondominant solution, then determine if the solution is not dominated by the solution in the external file Rep. If it is satisfied, add the solution to the external file Rep, and then delete the other solutions that are dominated by the solution.Third, IMBO algorithm was used to generate the next generation population. Repeat the above process until the stop condition is met.For the second step above, the process of storing the nondominant solution is described in detail below. Figures1 and 2 show the storage of nondominant solutions when the external file Rep is empty and not empty, respectively.Figure 1
Storage when the external file is empty.Figure 2
Storage when the external file is not empty.Figure1 shows the storage situation when the external file Rep is empty. At this point, it is directly judged whether the solution is a nondominant solution. If it is not, it is directly discarded (as shown in Case 1). If the solution is nondominant, it is added directly to the external file Rep (as shown in Case 2). Figure 2 shows the case when the external file is not empty. In this case, the first step is to determine whether the solution is a nondominant one. If the solution is not a nondominant solution, it is discarded directly (as shown in Case 1). If the solution is a nondominant solution, it needs to be further determined whether it is dominated by an existing solution in the external file Rep. If the solution is dominated by the solution in an external file, it is discarded (as shown in Case 1). If the solution is not controlled by an existing solution in the external file, it is added to the external file Rep. At this time, delete the other solutions in the external file that are dominated by this solution. Case 2 in Figure 2 above indicates that there is no other solution governed by this solution in the external file, so it is added directly to the external file Rep. Case 3 indicates that x2 and x3 are dominated by the new solution after adding the new solution to the external file. In this case, x2 and x3 need to be removed from the external file.When the MBO algorithm is updated iteratively, the corresponding dimension values of randomly selected individuals are directly assigned to the new individuals, which limits the diversity of the population and is detrimental to its evolution. For this reason, some ideas of differential evolution and particle swarm optimization (PSO) are added to the MBO algorithm to increase the diversity of the population. In differential evolution algorithms, the values on the corresponding dimensions of two individuals are subtracted to produce a difference vector, which is assigned a weight and added to another vector to produce a variation vector, thus increasing the diversity of the population. PSO algorithms update iteratively according to the individual extremes and global extremes, which can speed up convergence. Based on the above analysis, add the following actions to the migration operation:(18)xj,kt+1=xpbest,kt+xneighbor1,kt−xneighbor2,kt×R.In the above formula,xj,kt+1 is the k-th dimension of xj at the t+1 iteration, xpbest,kt is the k-th dimension of the individual extremum at the t-th iteration and xneighbor1t and xneighbor2t are two different individuals randomly selected from Subset2 during the t-th iteration. R is the weight assigned to the difference vector and can be calculated from the following formula:(19)R=rand−0.5×2.As can be seen from the above formula,R is a random number between −1 and 1. The larger the R, the greater the influence of the difference vector on the variation vector.NSGA-II uses a fast nondominant sorting method, which reduces computational complexity. Adding some operations in NSGA-II to the MBO algorithm can achieve better performance.During the iteration, the populations are merged after the migration and adjustment operations. The combined population was placed in a nondominated order, and the crowding degree of each individual was calculated. Then, the parent individuals are selected according to the crowding degree, and the solution is updated by the crossover operator. The following formula is the calculation method of new individualsXA′ and XB′ after the arithmetic crossover operation of XA and XB:(20)XA′=aXB+1−aXA,XB′=aXA+1−aXB.In the formula above,a is a random number between −0.5 and 1.5. The above operation can better maintain the diversity of the population.The flow chart of IMBO algorithm is shown in Figure3. The specific implementation process is as follows: first, initialize the population, external file Rep, population size, and other parameters. Then, the fitness values of all individuals in the population are calculated, and the nondominated solutions are found and stored in Rep; record the best position of the individual pbest, and the global optimal solution gbest. Finally, the algorithm is executed iteratively. The population is divided into two subpopulations, namely subset1 and subset2, which are transferred and adjusted to generate new solutions. Two new subpopulations are merged and nondominated sorting is performed. The crowding degree is calculated at the same time. Then, the parent individuals are selected for arithmetic crossover operation to further update the solution. The fitness of each individual is calculated and Rep, pbest, and gbest are updated. When the termination condition is satisfied, the algorithm is terminated; otherwise, continue the iteration.Figure 3
Schematic flowchart of IMBO.


## 4. Experimental Validation

Specifications of machining model are shown in Tables1 and 2.where Dw is the diameter of the workpiece, ap is the cutting depth, TR is the tool change time, Z is the number of milling cutter teeth, Nmin is the minimum spindle speed, Nmax is the maximum speed, η is the rotation coefficient of the machine tool, rε is the radius of tool tip arc; Ramax is the maximum surface roughness. All variables are explained in Chapter 2. According to the cutting manual, the wear resistance coefficient and cutting force coefficient of the cutting tool can be found, as shown in the table below:Table 1
Model specifications.
Dw/mmap/mmTR/sZrε/mmPcmax/kW·hFfmax/Nη10026004120080000.8RamaxL/mmNmin/r∙min−1Nmax/r∙min−1vfmax/m∙min−1vfmin/m∙min−1ae3.231594560008000360Table 2
Tool wear coefficient and cutting force coefficient.
Tool wear coefficientCutting force coefficientCvkhmhqhyhuhphCFqFxFwFyFKFCuF10670.10.20.250.20.10.1579001.31.00.20.750.251.1The material selected in the experiment is Ti-6Al-4V titanium alloy plate. According to [21], the cost parameter Ct is set to 60, the sum of labor cost C1 and management cost C0 per unit time is 2.39 yuan/min, the tool change time Tc of each process is 2.39 minutes, and the auxiliary time Tot is 0.1 minutes.The IMBO algorithm is simulated in MATLAB 2017b. In the algorithm, the population size of monarch butterfly is set to 120, and the number of iterations is set to 50. Compared with NSGA-II algorithm [4], PSO algorithm [5, 21] and traditional MBO with the same population size and iteration times, as shown in the figure below, IMBO algorithm has better optimization effect.As can be seen from Figure4, the IMBO algorithm has better optimization results than NSGA-II and traditional MBO algorithm. Here are the last six sets of data.Figure 4
Pareto diagrams of the four algorithms.According to Tables3 and 4, when the processing time is short, the production efficiency is high, but the tool wear is serious at this time, so the production cost is high. Extending the processing time can reduce the tool wear, but also reduce the production efficiency. When the processing time reaches 10.3 minutes, the production cost is not significantly reduced by extending the processing time.Table 3
Processing time/min.
IMBO10.3310.3510.3710.4110.4510.50NSGA-II10.2610.2810.3310.3610.4210.45MBO10.2910.3410.3510.3910.4510.47PSO10.2810.3310.3510.3810.4410.49Table 4
Cost/yuan.
IMBO26.2026.1526.1226.1026.0826.05NSGA-II27.1926.7926.4926.4026.3126.25MBO26.6926.5026.3026.2126.1526.14PSO26.6726.4326.3526.2326.1726.13The comprehensive evaluation method for remanufacturing of mechanical and electrical products is used to analyze six groups of data. The method is an analytic hierarchy process. The weight is determined by layer-by-layer refinement and level-by-level comparison, and then the weight of each subtarget to the total target is synthesized. The value of processing time function and production cost function is taken as the criterion layer, and the optimal parameter result is taken as the following objective layer:(21)T=11/41/41/51/21/34111/2324111/23252214321/31/31/411/231/21/21/321,C=111411/2112411/211/21531/21/41/41/511/31/3111/3311222311.The maximum eigenvalues are 6.07 and 6.35, respectively, and the weight vector of criterion layer to target layer isW2=0.5,0.5T. Consistency index CI=0. The eigenvectors corresponding to the maximum eigenvalues of T and C are obtained, and then the weight vector W1 from the scheme layer to the criterion layer is obtained by normalization as follows:(22)W1=0.0490.160.2060.190.2060.190.3420.050.0750.120.1220.30,W=W1×W2,=0.10450.19800.19800.19600.09750.2110.According to total orderW, the last group has the largest weight. That means when the cutting speed is 12.91 m/min, the feed rate is 0.16 mm/Z and the processing time of a single part is 10.41 min, the unit cost is 26.27 yuan. At this time, the production efficiency and production cost can reach a better value. It can be seen that different milling parameters have an important impact on production efficiency and production cost.


## 5. Conclusion

In this article, the mathematical optimization model of milling parameters is established with the goal of minimum production cost and maximum production efficiency, and the constraint violation degree method is used to deal with the constraints considering the constraints in actual processing. Some ideas of differential evolution algorithm and particle swarm optimization algorithm are added to the MBO algorithm, which can effectively speed up the convergence speed and avoid falling into the local optimal value. The IMBO algorithm is used to solve the optimization model and get a better parameter combination. Compared with NSGA-II and other algorithms, it shows that the IMBO algorithm can deal with the cutting parameter optimization problem well, and it also shows that IMBO algorithm has some advantages in dealing with multiobjective optimization problems [23].



---

*Source: 1348624-2023-01-11.xml*