Modified honey bees mating optimization algorithm for multi-objective uncertain integrated process planning and scheduling problem

The new technology of intelligent manufacturing makes the data of process planning and shop scheduling easier to interconnect, and the integration optimization of different manufacturing processes is an important technology to ensure the implementation of intelligent manufacturing. Integrated process planning and scheduling is a significant research focus in recent years, which could improve the performance of manufacturing system. At present, the research on integrated process planning and scheduling is insufficient to consider the multi-objective and uncertain characteristics widely existing in real manufacturing environment. Therefore, multi-objective integrated process planning and scheduling problem with uncertain processing time and due date is addressed in this article. The mathematical model of multi-objective uncertain integrated process planning and scheduling problem with uncertain processing time and fuzzy due date is established based on fuzzy set, in which the calculation method of uncertainty measurement objective is designed. An effective modified honey bees mating optimization algorithm has been designed to solve the proposed model. Queens set is constructed to maintain the non-dominated solutions found in the optimization process. The calculation method of mating probability between drone and queen bee based on Euclidean distance is designed. Fuzzy operators were utilized to evaluate fitness, judge the non-dominated relationship, and decode the scheduling solution. Different instances were designed and carried out to test the performance of the proposed method. The results show that the proposed method is very effective for solving multi-objective uncertain integrated process planning and scheduling.


Introduction
Modern manufacturing faces multiple pressures from the economy, society, and the environment, and manufacturing companies need to further improve efficiency, stabilize production, and enhance brand effects. Research on the theory of operational optimization of advanced manufacturing systems is particularly important for the development of manufacturing companies. In recent years, the rapid development of technologies such as intelligent manufacturing and industrial robot systems has greatly shortened the gap between theoretical research and practical applications. Process planning and shop scheduling are very important in modern manufacturing system, which have vital impact on the product processing capabilities, resource utilization, and production efficiency of the manufacturing system. The purpose of process planning is to select the appropriate processing methods and processing parameters for jobs by analyzing the machining characteristics, accuracy, materials, and so on. 1,2 After the process plans of jobs are determinate, the scheduling is to allocate the operations of all these jobs on machines over time by satisfying the precedence constraints in the process plans. 3 Although they are independent of each other and have different functions, they are mutually restricted. Considering process planning and shop scheduling as separate subsystems is not conducive to improve the production efficiency and equipment utilization ratio of the manufacturing system, which will bring a series of conflicts. 4 If only the operation sequence and processing resources of a single job are taken into account in the process planning stage, it is very easy to generate resource conflicts in the shop scheduling process. To match the real-time status of the workshop, 20%-30% of the previously processed plans must be replanned. 5 With the rise of intelligent manufacturing, the new technology of intelligent manufacturing makes the data of process planning and shop scheduling easier to interconnect, and the integration optimization of different manufacturing processes is an important technology to ensure the implementation of intelligent manufacturing. [6][7][8] Integrated process planning and scheduling (IPPS) can effectively overcome the conflicts caused by independent research, thus improving the efficiency of manufacturing systems. 9 IPPS problem is more complicated compared with process planning and shop scheduling problems which are non-deterministic polynomial (NP)-complete problems. 10 The IPPS problem has complex machine flexibilities, operation sequence flexibilities, and processing flexibilities, which makes problem solution space enormous and intricate. 11 As a result, developing an effective algorithm to deal with IPPS is a challenging work. Furthermore, many different objectives are involved in IPPS problem, such as makespan, total machine workload, total flow time, total tardiness, and so on. 12,13 Decision makers need to make a trade-off among different objectives. On the other hand, there are often a series of uncertain factors in actual production scheduling, such as uncertain processing time for operations, breakdown of machinery, due date variation, order change, and so on, 14 which makes it difficult for the deterministic optimization scheduling scheme to be executed smoothly during implementation. Srinivasan et al. 15 pointed out that the model mismatch caused by the existence of uncertain factors has become the bottleneck of application optimization technology in actual production. Considering multi-objectives and uncertain factors simultaneously could make IPPS problem more representative of real-world situations.
In recent years, researchers begin to study IPPS problem under uncertain environments, but only a few of these studies pay attention to multi-objective uncertain IPPS problem. A large number of scholars have used fuzzy set theory to study uncertain job shop scheduling problems, focusing on the uncertain processing time described by triangular fuzzy numbers and the uncertain due date described by triangle fuzzy numbers. 16 The research results of uncertain job shop scheduling problem can be extended to multi-objective uncertain IPPS problem. Swarm intelligence optimization algorithms based on Pareto optimization have been widely utilized to settle various multi-objective shop scheduling problems. [17][18][19][20] Honey bees mating optimization (HBMO) algorithm has been proposed based on the marriage behavior of the bees, which simulated the mating process of the queen of the hive. It has been applied to many combinatorial optimization problems successfully. 21,22 It is a meaningful work to apply HBMO to multi-objective uncertain IPPS problem based on fuzzy set and Pareto optimization.
In this article, multi-objective IPPS problem with uncertain processing time and due date is addressed by modified HBMO (MHBMO). The mathematical model of multiobjective uncertain IPPS problem with uncertain processing time and fuzzy due date is established based on fuzzy set, in which the calculation method of uncertainty measurement objective is designed. The following operations are designed to modify HBMO to solve multi-objective uncertain IPPS problem. Since the result of multiobjective optimization problem is a set of solutions rather than a solution, the queens set is constructed to maintain the non-dominated solutions found in the optimization process. Because of the multiple objectives that need to be considered, the calculation method of mating probability between drone and queen bee based on Euclidean distance is designed. The brood differentiation strategy is designed based on fast non-dominated sorting scheme. The instances of multi-objective uncertain IPPS problems with different scales were designed, and the proposed optimization method is utilized to solve the instances. The experimental results verify the effectiveness of the proposed method.
The remainder of this article is organized as follows. The "Literature review" section reviews the related works about IPPS in recent years. Problem description is elaborated in the "Problem formulation" section. The "Proposed MHBMO for multi-objective uncertain IPPS" section presents the workflow of the proposed MHBMO for multi-objective uncertain IPPS. The "Experiments and discussion" section shows the experiment results and discussion while the conclusion and future works are given in the last section.

Literature review
In manufacturing companies, process planning and shop scheduling are often two separate departments. It is not realistic to integrate the two departments. Considering this practical problem, IPPS mainly focuses on how to increase the information interaction between process planning and shop scheduling. According to different information interaction methods, the integration model of process planning and shop scheduling is mainly divided into three categories 1 : nonlinear method, closed-loop method, and distributed method. There are also some novel IPPS models proposed in recent years. Chu et al. 23 proposed a novel method for IPPS under production uncertainties. The integrated problem was formulated into a bi-level program. To solve the integrated problem, a hybrid method was developed, which iterated between a mixed-integer linear programming solver for the planning problem and an agent-based reactive scheduling method. Jin et al. 24 developed some novel mixed integer linear programming (MILP) models for IPPS in a job shop flexible manufacturing system, which were able to express and utilize flexibilities contained in network graphs, and hence had the power to solve network graph-based instances.
The established models had been tested on typical test bed instances to verify their correctness. Rietz et al. 25 addressed the integration of medium-term production planning and short-term scheduling models. A new detailed scheduling model based on a pseudo-polynomial network flow formulation that can be used to solve real size instances was described and analyzed. Different strategies were explored to simplify the model and reduce its number of constraints. Aguirre and Papageorgiou 26 presented a medium-term optimization-based approach for the integration of production planning, scheduling, and maintenance. A multiproduct single-stage batch process plant with parallel units and limited resources was considered in this research. Li 27 studied IPPS model for sustainable manufacturing and remanufacturing. Dai et al. 28 established an energy-saving IPPS model, in which the energy consumption objectives included the total energy consumption of all working machines in the manufacturing system during the production and nonproduction phases. Zhang et al. 29 proposed an integration model based on nonlinear process planning to implement energysaving method, and the Therblig-based model was used to predict the energy consumption of machine tools in product manufacturing process. It can be seen from the above studies that more optimization objectives and constraints were considered in the establishment of the IPPS model to provide effective serve for the actual production operation process.
According to the different number of optimization objectives considered during the search process, the IPPS solution method can be divided into single-objective optimization method and multi-objective optimization method. Li et al. 9 presented an effective hybrid algorithm based on genetic algorithm and variable neighborhood search to solve IPPS problem in packaging machinery factory. Zhang and Wong 30,31 proposed an object-coding genetic algorithm and constructive meta-heuristics to deal with IPPS problem. Zhang et al. 32 combined extended imperialist competitive algorithm with genetic algorithm to solve the distributed IPPS problem. Ausaf et al. 33 designed a priority-based optimization algorithm to settle multi-objective IPPS problem. Jin et al. 34 proposed a multi-objective memetic algorithm for IPPS problem. Zhang et al. 35 studied the IPPS problem in remanufacturing environment and proposed a multiobjective genetic algorithm based on simulation to optimize the optimization objectives in process planning stage and shop scheduling stage. Luo et al. 12 proposed a multiobjective genetic algorithm based on immune principle and external archive to settle multi-objective IPPS. Manupati et al. 36 constructed a mobile agent-based approach to optimize the IPPS problem with maximum completion time and machine utilization for the network manufacturing environment. Milica et al. 37 proposed a particle swarm optimization algorithm based on chaos theory to optimize the maximum completion time, mean flow time, and machine load balance of IPPS. Elahe 38 used the weighted method to deal with three objectives of the IPPS problem and solved the problem using a genetic algorithm.
In general, the research method of single-objective IPPS problem is more mature, mainly focusing on the design of effective hybrid algorithms. Multi-objective IPPS solution methods can be divided into non-Pareto methods and Pareto-based methods. The non-Pareto methods belong to the pretest method, that is, the determination of the weight represents the preference of the decision maker, and the setting of the weight will bring certain subjective influence factors to the optimization result. The Pareto-based methods belong to the posttest method, that is, the optimization is made first. The optimization process is only related to the characteristics of the problem itself, and less subjective factors. Because the Pareto-based method can more objectively reflect the characteristics of the problem, it is the focus of current research. This article will adopt Paretobased methods to handle multi-objective uncertain IPPS problem.
To make theoretical results better serve actual production, more and more researchers have begun to pay close attention to IPPS under uncertain environment. Seker et al. 39 applied genetic algorithm and fuzzy neural network to solve the uncertain IPPS problem in mass customization environment. The genetic algorithm was used to obtain the optimal scheduling sample data, and the fuzzy neural network was adopted to obtain the dynamic scheduling knowledge to guide the regeneration of the optimal scheduling scheme in the dynamic workshop environment. Liu et al. 40 proposed the ant colony optimization algorithm to solve the IPPS problem by considering the uncertainty of the arrival of new jobs. Xia et al. 41 constructed a mathematical model of dynamic IPPS problem, considering the dynamic events of machine faults and new job arrivals. The rolling window technology was applied to the solution of dynamic IPPS problems. Yu et al. 42 proposed a hybrid geneticparticle swarm optimization algorithm to solve the dynamic IPPS problem. The solution process was divided into static phase and dynamic phase. In the dynamic phase, IPPS could better cope with dynamic events by adjusting the process plans of jobs to be processed during the machine failure time period. Wen et al. 43 designed an improved genetic algorithm to settle multi-objective uncertain IPPS model by considering fuzzy processing time. Wang et al. 44 used interval numbers to describe the uncertainty of processing time in IPPS problem.
The research on IPPS problems under uncertain environment focuses on the following two aspects: using fuzzy numbers and interval numbers to describe the uncertainty of processing time and due date, and using rolling window technology, fuzzy neural network, increasing buffer time window, and so on to respond the uncertain events, such as machine failures and new job arrivals. Compared with IPPS problem in static environment, the uncertain IPPS problem brings the new uncertainty evaluation objectives. Due to the randomness of the uncertain events, the calculation of uncertainty evaluation objectives for uncertain IPPS is very complicated. Therefore, this article will design a novel calculation method of uncertainty evaluation objective.

Fundamental theory of fuzzy sets
In actual production process, factors such as the condition of different machines in the workshop and the different degree of operators' proficiency have certain uncertainties, so the processing time of jobs is often difficult to express with a certain value. The customer's due date requirement for jobs is usually an interval, and the customer will have the highest satisfaction when jobs can be delivered within this interval. Early or delayed delivery will reduce customer satisfaction. According to the characteristics of processing time and due date in the actual production process, many researchers use the triangular fuzzy numbers to describe the uncertain processing time of jobs and trapezoidal fuzzy numbers to describe the uncertain due date of jobs in the research of uncertain job shop scheduling. 45 The scheduling type of multi-objective uncertain IPPS problem studied in this article is job shop scheduling, therefore, the triangular fuzzy number is used to describe the indeterminate processing time of jobs, and trapezoidal fuzzy number is used to describe the indefinite due date of jobs in mathematical model.
When solving the multi-objective uncertain IPPS problem, the following fuzzy number operations need to be used.
(1) Addition operation When calculating the fuzzy processing time and the fuzzy completion time of each job, the addition operation of the triangular fuzzy numbers need to be used. If the two triangular fuzzy numbers areM ¼ ðm 1 ; m 2 ; m 3 Þ andÑ ¼ ðn 1 ; n 2 ; n 3 Þ, the addition operation is given in formula (1) (2) Taking large operation In the decoding process of the job shop scheduling stage, the fuzzy start time of each operation need to be computed by using the taking large operations of the triangular fuzzy numbers. If two triangular fuzzy numbers arẽ M ¼ ðm 1 ; m 2 ; m 3 Þ andÑ ¼ ðn 1 ; n 2 ; n 3 Þ, the taking large operation is given in formula (2) After the two triangular fuzzy numbers having undergone the taking large operation, the result is an approximate triangular fuzzy number.

(3) Comparison operation
Comparison operation is needed during decoding process of job shop scheduling, judging the dominant relationship and determining the fuzzy maximum completion time. If the two triangular fuzzy numbers areM ¼ ðm 1 ; m 2 ; m 3 Þ andÑ ¼ ðn 1 ; n 2 ; n 3 Þ, the following criteria are adopted to rank fuzzy numbers: and if C 3 ðMÞ > C 3 ðÑ Þ, thenM >Ñ .

Mathematical model for multi-objective uncertain IPPS
The multi-objective uncertain IPPS problem focused in this article can be defined as follows 11 : There are a set of n jobs to be processed on m machines. Each job has various operations and alternative process plans. The processing time and due date of each job is uncertain. One operation is allowed to be executed by alternative machines. The aim of multi-objective uncertain IPPS is to select suitable manufacturing resources for each job, determine the operations' processing sequence and the start time of each operation on each machine by satisfying the precedence constraints among operations and achieving the corresponding multiobjectives simultaneously.
Based on the multi-objective IPPS mathematical model 12 and the basic theory of the above fuzzy sets, a multi-objective uncertain IPPS optimization problem model is established. To better describe the model, the following symbols are defined at first: n: the number of all jobs to be machined; m: the number of all machines in the production workshop; D i : the due date of job i, which is descried as trapezoidal fuzzy numbers; g i : the number of all alternative process plans for each job; n il : the number of operations contained in the lth process plan of job i; o ijl : the jth operation in lth process plan of job i; k: the optional machine number of o ijl ; P ijlk : processing time of o ijl machined in machine k, triangular fuzzy numbers; C ijlk : the earliest completion time of o ijl machined in machine k, triangular fuzzy numbers; C i : the completion time of job i, triangular fuzzy numbers; A: a very large positive number; The optimization objectives involved in multi-objective uncertain IPPS problem can be divided into two categories. One is the objectives represented as triangular fuzzy numbers, such as fuzzy maximum makespan, fuzzy maximum machine load, and fuzzy total machine load. Another is the objectives about the uncertainty measure, which contains the mean customer satisfaction and the spread of makespan. The five objective functions are provided as follows: f 2 : Minimize fuzzy maximum machine load f 3 : Minimize fuzzy total machine load f 4 : Maximize average customer satisfaction ST i ð7Þ The following constraints need to be satisfied: (1) The earliest completion time of the first operation in the lth process route of job i need to satisfy the following constraint (2) The earliest completion time of the last operation in the lth process route of job i need to satisfy the following constraint (3) Operation constraint: Only one operation of each job can be processed at the same time (4) Machine constraint: Each machine can only process one process at a time i;x2½1;n;j;y2½1;n il;xz ;l;z2½1;g i;x ;k2½1;m ð12Þ (5) Process plan constraint: Only one process plan can be selected each job at a time (6) Only one optional machine can be selected for each operation X m k¼1 Z ijlk ¼ 1 i 2 ½1; n; j 2 ½1; n il ; l 2 ½1; g i ð14Þ (7) The completion time of all operations should meet the following constraint In the above multi-objective uncertain IPPS problem model, the objectives for uncertainty measure contain mean customer satisfaction and the spread of fuzzy makespan. The greater the average customer satisfaction, the more satisfied the customer is with the uncertain due date. The smaller the spread of the fuzzy maximum completion time, the smaller the uncertainty of the fuzzy makespan. The evaluation method of customer satisfaction in multiobjective uncertain IPPS problem will be discussed in the following section.

Evaluation method of customer satisfaction in uncertain IPPS
In a deterministic environment, makespan and due date of the job are determined value, and it is easy to evaluate whether the job is delayed or completed in advance. However, when the makespan and due date of the job are uncertain values, the relationship between uncertain makespan and the due date of each job is very complicated. How to evaluate the relationship between the makespan and the due date of the job is very important. In the research of uncertain job shop scheduling problems, scholars generally use customer satisfaction to evaluate the relationship between fuzzy makespan and fuzzy due date. It is very necessary to consider optimizing customer satisfaction when conducting a research on multi-objective uncertain IPPS problem.
Suppose in a scheduling scheme, the fuzzy makespan of andD i , respectively, are shown in formulas (16) and (17) Traditionally, most of the literature used customer satisfaction calculation method based on the intersection area of membership function, which is given in formula (18) The method has the following problems: (1) According to the different positions between the three vertices of the triangular fuzzy number and the four vertices of the trapezoidal fuzzy number, relationships will be generated between the fuzzy completion time and the fuzzy delivery date. The calculation method of the intersection area of each position relationship is generally different, which will greatly increase the calculation time of the customer satisfaction during the problem-solving process.
(2) When the makespan of job i is located at three different positions as shown in Figure 1, according to the traditional customer satisfaction evaluation method, customer satisfaction is the same in all three locations. However, as can be clearly seen from Figure 1, when the makespan of the job i is in the position 2, all possible makespan can ensure maximum customer satisfaction of the job. When the completion time of job i is at position 1, if the makespan of job i is between (d i1 , d i2 ), the customers' satisfaction cannot be maximized. Similarly, when the completion time of job i is at the position 3, the makespan of job i may drop in (d i3 , d i4 ), and the customer satisfaction cannot be maximized at this time. It can be seen that the decision makers tend to choose the scheme in the actual production Figure 1. Three different positional relationships between fuzzy makespan and fuzzy due date.
process when the makespan of the job i is at the position 2. Compared with the other two schemes, position 2 can avoid the risks in the uncertain environment. However, if formula (18) is used as the customer satisfaction evaluation index, it is impossible to distinguish these three different schemes.
If the decision maker evaluates the customer satisfaction based on formula (18), it cannot evade the deterioration of the objective value in the worst case of the uncertain environment.
In view of the shortcomings of the evaluation method based on formula (18), this article makes some improvements on the satisfaction evaluation method proposed by Hu et al., 46 which is used to comprehensively evaluate the customers' satisfaction with the scheduling scheme under the uncertain environment. The proposed method is consistent with the evaluation method based on intersection area, while the calculation is relatively simple. This method of evaluating customer satisfaction is detailed as below.
Since the makespan of job i is a triangular fuzzy number, it is necessary to consider not only the degree of customer satisfaction at the most possible completion time but also the influence of customer satisfaction when the completion time deviates from the most likely completion time. As shown in Figure 2, D i c i2 ð Þcould be directly used to evaluate customer satisfaction with the most possible completion time of job i. However, when evaluating the customers' satisfaction with the possible completion time of the job i in the worst case, it is not suitable to directly use D i c i1 ð Þ and D i c i3 ð Þ for evaluation. This is because c i1 and c i3 represent the most optimistic completion time and the most pessimistic completion time of job i, that is, the worst possible result, but their credibility is 0, so there is no reference value. Suppose the customer satisfaction is ST il when the worst possible result occurs while only the delay time of job i is considered. Suppose the customer satisfaction is ST ie when the worst possible result occurs while only earlier completion of job i is considered. In this case, the following definitions are given According to the possibility theory proposed by Zadeh, 47 ðÀ1;D i is less than or equal to the fuzzy number ofD i , ½D i ; þ1Þ is more than or equal to the fuzzy number ofD i . As shown in Figure 2 . According to the mathematical meaning, ST il represents the minimum possibility that the fuzzy completion time is less than the fuzzy due date. The greater the possibility that the fuzzy completion time is smaller than the fuzzy delivery date, the larger ST il is. As a result, ST il can be used to measure the customers' satisfaction with the completion of job i only the tardiness penalty is considered. Similarly, ST ie represents the minimum possibility that the fuzzy completion time is more than the fuzzy due date. The greater the possibility that the fuzzy completion time is larger than the fuzzy due date, the larger the ST ie is, so ST ie can be used to evaluate the customers' satisfaction with completion of job i only the early penalty is considered. As can be seen from Figure  credibility. They are able to assess the changes in customer satisfaction in the worst case while taking into account certain credibility.
In the case where only the job is considered to be penalized in advance, when the completion time of the job is greater than or equal to d i2 , the job may not be prematurely completed, and the customers' satisfaction is maximum. When the completion time of the job is less than d i2 , there is a possibility of early completion, and the customers' satisfaction is gradually reduced. In the case where only the job is considered to be penalized in tardiness, when the completion time of the job is less than or equal to d i3 , the job is unlikely to be completed, and the customers' satisfaction is maximum. When the completion time of the job is larger than d i3 , there is a possibility of tardiness, and the customers' satisfaction is gradually reduced.
Based on the above analysis, it can be concluded that there are three possibilities for the value of ST il : There are also three possibilities for the value of ST ie : The early or delayed completion of jobs will result in a certain degree of decrease in customer satisfaction. Under normal circumstances, the customers' tolerance for the early completion of the job is higher than that of tardiness. Therefore, in this research, ST il and ST il are given different weights ! l and ! e , while ! l > ! e ; ! l þ ! e ¼ 1.
Based on the above discussion, the customers' satisfaction with job i can be computed by formula (21) This evaluation method only needs to pay attention to the position of the three points (c ie , . This calculation method is simpler than the evaluation method of the intersection area, which can better avoid the impact of worst possible results on the solution in an uncertain environment and reflect the customers' satisfaction with a given solution.

Proposed MHBMO for multi-objective uncertain IPPS
The basic principle of HBMO The HBMO algorithm simulates the mating process of the queen in the hive. At the beginning of the whole mating flight, the queen flies away from the nest and performs a dance, and then the drones follow the queen and mate with her in the air. It means that some sperm from the drones will reach the queen's spermatheca to form the genetic pool of the colony. The queen uses the stored sperms to fertilize eggs. Each time a queen lays fertilized eggs, she retrieves a mixture of the sperms stored in the spermatheca randomly. The queen mates multiple times, and the drones will die after mating with the queen once. The mating flight can be considered as a set of transitions in a state-space. At the start of the flight, the queen is initialized with some energy and speed, and she will return to her nest when the energy is less than a threshold or when her spermatheca is full. The queen moves between the different states in the space in some speed and mates with the followed drones at each state probabilistically. A drone mates with the queen probabilistically using an annealing function as follows 21 where Prob(Q, D) is the probability of adding the sperm of drone D to spermatheca of queen Q (i.e. the probability of a successful mating); Dðf Þ is the absolute difference between the fitness of D and the fitness of Q; and Speed(t) is the speed of the queen at time t. It is apparent that the probability of mating is high either when the queen is still in the start of her mating flight or when the fitness of the drone is as good as the queen's. If the mating is successful, the drone's sperm is stored in the queen's spermatheca.
After each transition in space, the queen's speed and energy are reduced using the following equations 21 where the a 2 ð0; 1Þ, it means that the speed and energy of the queen will be reduced after each transition and each step.
After the mating flight, the queen uses the sperms stored in the spermatheca to create the new broods by the crossover operator. Because the broods contain different parts of the genotype in the drones' sperm, the quality of the new broods has more exploration abilities.
As the role of workers in the colony is restricted to brood care, the workers are employed to improve the broods produced by the crossover operator. In HBMO algorithm, each worker is represented as a heuristic which help improve a set of broods. The number of workers contains two parts (w ¼ w1 þ w2), w1 is the number of single local search strategies and w2 is the number of their combinations. Each new produced broods will choose one worker to improve it randomly. The queen will be replaced by the new brood, which has better fitness than her. Otherwise, the brood becomes one of the drones in the next flight if its fitness is better than one of the drones.
With the basic characters of the honey bees' mating described above, it is clear that the HBMO algorithm consists of three main components: 1. The queen's mating flight: In this component, the queen starts a mating flight and mates with the drones followed. 2. Generation of broods: In this component, a queen generates a number of new broods using the sperms in its spermatheca. 3. Improvement of broods: In this component, the new produced broods choose a worker to improve itself randomly. This component can increase the exploration abilities of the algorithm.
The basic HBMO algorithm is well applied in combinatorial optimization problems and is generally superior to genetic algorithm and particle swarm optimization algorithm. But it only works for single-objective optimization problems. To solve the problem that the basic HBMO algorithm is difficult to optimize multi-objectives at the same time, this article mainly improves the basic HBMO algorithm as follows: 1. There is only one queen in the basic HBMO algorithm, which represents the optimal solution found in the process of optimization. Since the target of multi-objective optimization problem is to find a set of non-dominated solutions rather than a solution, a set of queens is designed in MHBMO to save the non-dominated solutions found in the algorithm optimization process. 2. Each bee has multiple evaluation objectives when solving multi-objective optimization problems, so the calculation method of Dðf Þ in formula (22) needs to be changed. Therefore, a method for calculating the mating probability of drones and queen bees based on Euclidean distance is designed. 3. In the single-objective optimization problem, the differentiation pattern of the broods can be carried out by simply comparing the fitness value. This method is not suitable for multi-objective optimization problems. Aiming at the multi-objective optimization problem that needs to be solved in this article, a broods differentiation strategy based on fast non-dominated sorting is designed.
The details of the above specific changes to HBMO are described in specific in the "MHBMO for uncertain job shop scheduling" section.

The framework of the proposed method
Based on the integration model proposed by Li et al., 9 the workflow of the proposed MHBMO for multi-objective uncertain IPPS is given in Figure 3.
The main steps of the proposed method are described as follows: Step 1: Generate the initial process plan population for each job randomly.
Step 2: Optimize the process plan for each job by HBMO respectively.
Step 3: Output near-optimal process plans represented by queen for each job.
Step 4: According to the determinate process plan for each job, generate the initial job shop scheduling population randomly.
Step 6: Output the optimized process plan for each job and the corresponding scheduling plan.
Step 7: Update the non-dominated set for uncertain IPPS.
Step 8: If the number of current generation reaches to the maximum generation of uncertain IPPS, output the non-dominated set. Otherwise, go to step 2.

HBMO for uncertain process planning
The aim of process planning in this research is to provide various near-optimal process plans for the scheduling system. HBMO is applied to optimize process planning population in this article. Encoding and decoding schemes in Li et al. 9 were used to initial the population of process planning. The algorithm workflow is the same with the algorithm presented in Wen et al. 21 The crossover operations in Luo et al. 12 were utilized to generate broods. The local strategies as workers are the same with the method proposed in Wen et al. 21 As the processing times of operations are fuzzy numbers, the related operations of fuzzy numbers were utilized to evaluate fitness of bees in the population. The difference in fitness between different bees is calculated using formula (3).

MHBMO for uncertain job shop scheduling
After the uncertain process planning, fuzzy maximum machine load and fuzzy total machine load are determinate. The target of scheduling optimization in uncertain IPPS is minimizing fuzzy makespan, spread of fuzzy makespan, and customer satisfaction simultaneously. MHBMO is designed to tackle the multi-objective uncertain job shop scheduling problem. The flowchart of MHBMO for uncertain job shop scheduling is shown in Figure 4 and the details of the algorithm will be described as follows.
Initial population and fitness evaluation. For each individual in the scheduling population, the operation-based encoding method is used as the encoding method. 48 Assuming that the number of jobs to be machined is 4, after using the HBMO algorithm to optimize the flexible process plans of each job, the operation number of each job can be determined. Suppose that job 1 has six operations, job 2 has five operations, job 3 has three operations, and job 4 has five operations. A feasible encoding scheme in the shop scheduling population is [1 1 2 3 4 4 2 4 1 1 2 4 2 1 1 4 2 3 3]. The sixth element in the sequence is 4, 4 has been repeated twice, so this element represents the second operation of job 4. Each chromosome should be decoded into active schedules in the decoding procedure. 49 Since the operation-based encoding method does not generate infeasible solutions, this article randomly generates the initial bee colony of the job shop scheduling. The fuzzy makespan, spread of fuzzy makespan, and customer satisfaction of the scheduling scheme represented by each honey bee in the bee colony can be calculated after the decoding operation.
The fast non-dominated sorting approach in non-dominated sorted genetic algorithm-II (NSGA-II) proposed by Deb et al. 50 is utilized to divide the population into mutually disjoint non-dominated subpopulations (F 1 , F 2 , F 3 , . . . ) according to fuzzy makespan, spread of fuzzy makespan, and customer satisfaction. The individuals in the subpopulations with a high non-dominated rank are dominated by the individuals with a low non-dominated rank. The honey bees with rank 1 are the non-dominated individuals in the current population. After fast non-dominated sorting, all bees in the bee colony are given a non-dominated grade, and individuals with lower nondominated grades have better fitness values. Improvement of broods by workers. In the proposed MHBMO algorithm, each worker bee is equivalent to a local search strategy. Therefore, worker take care of broods is equivalent to using a local search strategy to update the broods. Since the job shop scheduling considered in uncertain IPPS is a multi-objective optimization problem, the calculation time required for each comparison between individuals is greatly increased compared to the single-objective optimization problem. Therefore, it is not suitable to spend too much time in the local search process. Therefore, a worker was randomly selected to improve broods.
Based on the encoding method used in the job shop scheduling stage, two neighborhood structures (N1, N2) are used to construct the workers, and the number of workers is two. The specific neighborhood structures are as follows: Based on the above strategy, workers can be compared to the local search process. The detailed steps are as follows: Step 1: Set the maximum number of iterations of workers to improve broods in job shop scheduling system is LS max , set t ¼ 1.
Step 2: Randomly select a worker from the candidate worker bee w i , i 2 f1; 2g. w i is employed to improve brood t, new brood t 0 will be obtained after the improvement.
Step 3: If brood t 0 can dominate brood t, t is replaced by t 0 . Step 4: Set t ¼ t þ 1, if t < LS max , go to step 2.
Otherwise, the cultivation operation of workers is terminated.
MHBMO for uncertain job shop scheduling. The main steps of the proposed MHBMO for uncertain job shop scheduling are as follows: Step 1: Initialize the parameters of MHBMO algorithm, including the number of iterations Gen max , attenuation coefficient of queen's velocity and energy a, energy threshold of the queen's mating flight threshold, the capacity of queen's spermatheca SperSize, number of honey bees PopSize, number of queens QueenSize, number of broods BroodSize, number of workers WorkerSize, and number of iterations of workers L max .
Step 2: According to the output of the process planning system, randomly generate PopSize honey bees according to the method in the "Initial population and fitness evaluation" section. Each honey bee represents a scheduling scheme. Calculate the three objective function values represented by each honey bee. Clear the queens set, set Gen ¼ 0.
Step 3: Initialization of drones population and queens set.
Step 3.1: Deposit non-dominated individuals in the current population into the queens set. If the number of non-dominated solutions in the current population is greater than the number of queens, then calculate the crowded distance of each individual in the non-dominated solution set and put the individuals with large crowding distance into the queens set. The calculation method of individual crowding distance can refer to the literature. 50 Step 3.2: Remove the QueenSize individuals deposited in the queen of bees from the colony. Randomly generate QueenSize individuals and add these individuals to the bee colony. At this time, PopSize individuals are used as drones population.
Step 4: Queen's mating flight. Randomly select a bee from the queens set as the queen, initialize the speed and energy of the queen, and empty the queen's spermatheca. Set a represent the number of drone genotypes in the queen's spermatheca, set a ¼ 0. Let t represent the number of flights of the queen, set t ¼ 0. If the energy of queen Energy ðtÞ > threshold and a < SperSize, repeat the following steps: Step 4.1: Randomly select a drone from the population.
Step 4.2: Calculate the probability value of the drone mating with the queen bee by formula (22).
Because the first objective value for each individual in the current population is triangular fuzzy numbers, the second objective and the third objective values are real numbers, Dðf Þ is calculated by formula (25) where Q f i and D f i represent the ith objective function value of the scheduling scheme represented by the queen bee and the drone, respectively. Randomly generate a real number Use the other individuals in P temp to update the drones population between 0 and 1, if the real number is less than the calculation result of formula (22), deposit the genotype of the drone into queen's spermatheca. Remove the drone from the drone colony. Set a ¼ a þ 1, otherwise, no operation is performed.
Step 4.3: Update the speed and energy of the queen according to formulas (23) and (24). Set t ¼ t þ 1.
Step 5: Broods generation. Randomly select a genotype from the queen's spermatheca, and use the crossover operators in the literature 21 to cross the queen to generate broods until BroodSize broods are generated.
Step 6: Improve the genotype of each brood using the operations designed in the "Improvement of broods by workers" section.
Step 7: Update queens set and drones population.
Step 7.1: Combine the original drones population and the broods population into a population P temp , divide P temp into non-dominated subpopulation F ¼ fF1, . Otherwise, return to step 7.4.
Step 8: The drone population is updated using the diversity maintenance strategy based on the immune principle proposed in Luo et al. 12 Step 9: Gen ¼ Gen þ 1. Determine whether the algorithm termination criterion is satisfied, and the algorithm termination criterion is to reach the maximum number of iterations of the algorithm. If Gen < Gen max , jump to step 4. Otherwise, output the final set of queens obtained from the scheduled population, each of which contains the process plans of each job and the final scheduling scheme.

Experiments and discussion
The above MHBMO algorithm for solving multi-objective uncertain IPPS problems is programmed by Visual Cþþ. The computer performance of the running program is 2.0 GHz Core (TM) 2 Duo CPU, 2.00 GB memory. Two sets of multi-objective uncertain IPPS instances are designed in this article. The proposed MHBMO was utilized to solve the instances. The value l of the diversity maintenance strategy based on the immune principle is set to 0.9, the non-dominated solution set size is set to 40, and the maximum number of iterations of the algorithm is set to 100. The algorithm termination criterion is that the current iterations number reaches the maximum iterations number of the algorithm. The algorithm parameters of other process planning and shop scheduling sections are shown in Table 1.

Experiment 1
Experiment 1 contains two instances (5 Â 5, 10 Â 5) of different scales. Among them, the processing information of the instance 5 Â 5 (five jobs processed in five machines) is shown in Tables 2 and 3. The first five jobs and the last five jobs of instance 10 Â 5 are the same with jobs in instance 5. To verify the performance of the proposed MHBMO algorithm for solving multi-objective uncertain IPPS problems, NSGA-II-based optimization method was used to solve multi-objective uncertain IPPS problem in experiment 1 simultaneously, the optimization method is the same as the process described in the fourth section, but the shop scheduling part is optimized using the NSGA-II algorithm. The two algorithms run independently 20 times, and the performance of the two algorithms is evaluated using generation distance (GD) and inverse GD (IGD). The calculation results are shown in Table 4. The relationship between the fuzzy makespan and the fuzzy due date of the five jobs with the mean satisfaction degree of 0.475705 and 0.288272 is shown in Figures 5 and 6, respectively.

Experiment 2
Experiment 2 contains three instances (instance 8 Â 8, 10 Â 10, and 15 Â 10) of different scales. In the instance 8 Â 8, the first five jobs have the same machining process and flexible process plans as five jobs of the instance 5 Â 5 in experiment 1. The last three jobs have the same processing information as the last three jobs in the instance 5 Â 5. Each operation has eight optional machines, in which the optional machine 1 to machine 5 of each process has the same fuzzy processing time as the five optional machines of instance 5 Â 5 in experiment 1. The fuzzy processing time of the optional machine 6 to the optional machine 8 is the same as the fuzzy processing time of the optional machine 3 to machine 5 in instance 5 Â 5. The fuzzy due date of the eight jobs is the same as the due date of the first eight jobs in instance 10 Â 5. The 10 jobs in instance 10 Â   10 have the same processing information as jobs in instance 10 Â 5. Each job has 10 optional machines. The machine 1 to machine 5 has the same fuzzy processing time as the machine 1 to machine 5 of the jobs corresponding to instance 10 Â 5. The fuzzy processing time of machine 6 to machine 10 is the same as the fuzzy processing time of machine 1 to machine 5. The fuzzy due date of the 10 jobs is the same as the fuzzy due date of the job in instance10 Â 5. The first 10 jobs in instance 15 Â 10 have the same machining information as the 10 jobs in instance 10 Â 10. The last five jobs have the same machining information as the last five jobs in instance 10 Â 10.The fuzzy delivery date of the first 10 jobs is the same as the fuzzy due date of 10 jobs in instance 10 Â 10, and the fuzzy due date of the last five jobs is the same as the fuzzy due date of the last five jobs in the instance 10 Â 10.
The MHBMO algorithm and the NSGA-II algorithm were used to solve the three different scale instances in experiment 2. Each algorithm runs 20 times independently. The performance of the two algorithms is evaluated by GD and IGD. The computation results are shown in Table 5.

Discussion
The computation results of the first example show that MHBMO can find more non-dominated solutions than NSGA-II in 20 independent running. At the same time, in the GD and IGD performance indicators, MHBMO is smaller than the NSGA-II in the maximum, minimum, and mean values. It can be seen intuitively from Figures 5 and 6 that when the average satisfaction is larger, there is a large intersection area between the fuzzy makespan of jobs and the fuzzy due date. It is shown that the customer satisfaction evaluation method designed in this article can guarantee the consistent evaluation results with the customer satisfaction evaluation method based on intersection area. The calculation results of experiments show that MHBMO can find more non-dominated solutions than NSGA-II in 20 independent running. Most of the maximum, minimum, and mean values of GD obtained by MHBMO algorithm are smaller than NSGA-II. All of the maximum, and mean values of IGD obtained by MHBMO algorithm are smaller than NSGA-II. In summary, MHBMO proposed in this article revealed better performance than NSGA-II in solving multi-objective uncertain IPPS problems. This   6.23841 Â 10 À3 Mean 6.66037 Â 10 À3 7.8177 Â 10 À3 100 non-dominated solutions were found, in which 58 were obtained by HBMO and 42 by NSGA-II. HBMO: honey bees mating optimization; MHBMO: modified honey bees mating optimization; GD: generation distance; IGD: inverse generation distance. The bold-faced values in Table 5 indicate that the values obtained by MHBMO are better than the values obtained by NSGA-II.
proposed method can effectively solve multi-objective uncertain IPPS.

Conclusion and future works
This article has conducted in-depth research on the multiobjective uncertain IPPS problem. Based on the fuzzy set theory, a multi-objective uncertain IPPS problem model is established, which considers uncertain processing time and uncertain due date. The objectives of the fuzzy objectives and the uncertainty of the metric scheduling scheme are taken into account comprehensively. The related fuzzy number operations were used to realize the fitness evaluation, non-dominated relationship judgment, and scheduling decoding. As the lack of multi-objective uncertain IPPS problem instances, instances of different sizes for multiobjective uncertain IPPS problem were designed. The designed instances were solved using the optimization method based on MHBMO and the optimization method based on NSGA-II. The experimental results show that the proposed method in this article can effectively solve the multi-objective uncertain IPPS problem and demonstrate better performance than NSGA-II.
With the development of fuzzy mathematics, some new fuzzy numbers are gradually proposed. It is an important research direction to consider the use of new fuzzy numbers to more accurately describe the uncertainty in IPPS. On the other hand, some novel multi-objective optimization algorithms with better non-dominated solution set construction method and diversity maintenance strategy can be studies to solve multi-objective uncertain IPPS.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.