Many-Objective Optimization and Decision-Making Method for Selective Assembly of Complex Mechanical Products Based on Improved NSGA-III and VIKOR

In Industry 4.0, data are sensed and merged to drive intelligent systems. This research focuses on the optimization of selective assembly of complex mechanical products (CMPs) under intelligent system environment conditions. For the batch assembly of CMPs, it is difficult to obtain the best combinations of components from combinations for simultaneous optimization of success rate and multiple assembly quality. Hence, the Taguchi quality loss function was used to quantitatively evaluate each assembly quality and the assembly success rate is combined to establish a manyobjective optimization model. The crossover and mutation operators were improved to enhance the ability of NSGA-III to obtain high-quality solution set and jump out of a local optimal solution, and the Pareto optimal solution set was obtained accordingly. Finally, considering the production mode of Human–Machine Intelligent System interaction, the optimal compromise solution is obtained by using fuzzy theory, entropy theory and the VIKOR method. The results show that this work has obvious advantages in improving the quality of batch selective assembly of CMPs and assembly success rate and gives a sorting selection strategy for non-dominated selective assembly schemes while taking into account the group benefit and individual regret.


Introduction
In the manufacturing environment of Industry 4.0, it becomes easy to accurately collect various data during production. Relying on a large amount of data, intelligent systems can optimize and make decisions in the production process. Most mechanical products are manufactured by the processing and assembly of components [1]. Under the condition of high requirements for product matching accuracy, due to the limitation of component processing capacity and manufacturing cost, it is unrealistic and uneconomical to completely rely on improving the accuracy of the processing process to meet and improve the product matching accuracy [2]. In practice, the demand of producers is to reduce manufacturing costs as much as possible while ensuring quality to obtain maximum product competitiveness. Therefore, selective assembly is one of the feasible methods to achieve high-precision assembly and reduce costs by using data and intelligent systems.
Researchers have made a lot of progress in selective assembly by researching the grouping of assembly components [3,4]. However, due to the many limitations of these methods, researchers have conducted further research on selective assembly methods from other perspectives. With the development of research, some intelligent optimization algorithms with better optimization effects have been used in the selective assembly field. Kannan et al. [5] proposed a genetic algorithm for selective assembly. In this research, the selective assembly process was analyzed, and a better combination was selected through the use of genetic algorithms to minimize assembly differences. Kannan et al. [6] used a genetic algorithm to obtain the best combination of the selection group with the smallest gap within the gap specification and at the same time minimized the deviation from the objective value by minimizing the Taguchi loss function. Fei et al. [7] proposed a new grouping method and chromosome structure using genetic algorithms to minimize the remaining components. Tan and Wu [8] summarized the research results of selective assembly and expanded the concept of selective assembly, distinguishing between direct selective assembly, using information directly obtained from measurements on component characteristics, and fixed bin selective assembly using components sorted into bins. The branch-demarcation algorithm is used to solve the selective assembly problem. Lu and Fei [9] proposed a new grouping method using a genetic algorithm to minimize the remaining components of a multidimensional chain and improved the genetic algorithm chromosome structure into a two-dimensional structure, making it more competitive in solving multidimensional chain problems. Manickam and De [10] proposed an optimal size selective assembly method, which is suitable for small batch quantities with wide tolerances, and successfully applied genetic algorithms in some components of the launch vehicle. Liu and Liu [11] proposed a method to determine the number of assembly groups for remanufactured engines. By improving the previous assembly method of specifically determining the number of groups, the number and interval of the components to be assembled are variable to reduce the waste of components. Wang [12] proposed the maximum assembly acceptance method under the condition of a non-normal distribution in matching components and optimized it by a genetic algorithm based on numerical experiments under different conditions. Kannan et al. [13] proposed a new selective assembly model to obtain specified tolerances in high-precision assembly and achieved a high success rate using genetic algorithms. Chu et al. [14] proposed an improved genetic algorithm to conduct application research on the selective assembly process of rotor vector (RV) reducers and verified the feasibility and effectiveness of the method. Jeevanantham et al. [15] studied the accumulation and extension of geometric and dimensional tolerances (GDTs) using tolerance analysis. By integrating GDT, the size deviation of components is grouped, and the optimal combination of component selective assembly is obtained by a genetic algorithm to obtain the minimum assembly deviation.
Kannan et al. [16] used the particle swarm optimization algorithm (PSO) to obtain the selective grouping and combination of component assembly, which significantly reduced the assembly deviation. Raj et al. [17] used the PSO algorithm to develop a code suitable for solving selective assembly problems and assembled components in batches to improve assembly efficiency. Rajesh et al. [18] studied the grouping selective assembly process, established the assembly loss model by using the Taguchi quality loss function, and obtained the combination of assembly losses that were as small as possible by using an improved immune algorithm. Rajesh et al. [19] proposed the Taguchi quality loss function based on symmetric intervals, studied the grouping selective assembly process, and used the improved sheep flock heredity algorithm to minimize the assembly clearance variation. Asha et al. [20] used a non-dominated sorting genetic algorithm (NSGA) to obtain the best combination of the selection group with the smallest change in the assembly clearance between the piston and the cylinder under the consideration of multiple characteristics. Although the multiple-characteristics-related situation is not considered in the assembly optimization process, it is of guiding significance for the selection of the processing and manufacturing process of the components. Raj et al. [21] applied the NSGA-II to batch selective assembly and applied it to complex components composed of pistons, piston rings, and cylinders. They proved that the algorithm is superior to existing methods in the literature in generating solutions with no remaining components and minimal gap changes. Researchers have applied NSGA in selective assembly and have begun to pay attention to solving the many-objective optimization problem in selective assembly.
Some of the latest research work has led to further development of selective assembly. Sun [22] applied the Taguchi theory to group and match the three-dimensional uniform transformation deviation, and the overall assembly quality was improved. Xing et al. [23], considering that the geometric error of the mating surface is an important factor affecting assembly quality, proposed an optimization method of shaft hole selective assembly based on relative entropy and dynamic programming, providing a new idea for the selective assembly of multiple batches of precision components. Rezaei Aderiani et al. [24] used digital twin technology to apply selective assembly methods to the assembly of automotive sheet metal components, which significantly improved the geometric variation and average deviation of the components. This research clearly discusses the difference between sheet metal assembly and linear assembly and notes that the selective assembly problem should be a many-objective optimization problem. Rezaei Aderiani et al. [25] improved the genetic coding method of intelligent optimization algorithms for solving selective assembly problems and proved in three cases that the effect of intelligent algorithms was better than mixed integer linear programming and mixed integer nonlinear programming in selective assembly, and this method could effectively reduce the calculation time. A methodology has been proposed in [26,27], which defines system requirements formally and automatically verifies and evaluates the system through simulation. Based on Formal Requirements Modeling Language and the Modelica Language, it is verified in two practical engineering project cases. This provides a theoretical basis and feasibility for formal confirmation of system requirements and simulation verification.
In engineering practice, the assembly of CMPs involves a multidimensional assembly dimension chain formed by multiple component sizes. From the literature survey, most of the above works focus on the improvement of selective assembly scheme of single dimension chain, less work has been done to model and research the quality loss of CMPs with a multidimensional assembly dimension chain, and there are certain limitations: (1) In most of works, the optimization objective in the selective assembly process is relatively single. The assembly quality loss simply sums up the quality loss of each separate assembly dimension chain without considering the independent quality of each dimension chain; (2) There have been very few previous studies considering the production mode of Human-Machine Intelligent System interaction and could not order and make decisions on the non-dominated solution set.
Therefore, this study aims to solve the problem of selective assembly in the manufacturing process of complex mechanical products, including quality optimization, cost reduction, and scheme decision making.
The proposed problem can be solved in three stages. Firstly, the selective assembly process of CMPs is described, and then this mathematical model is established using the success rate and the Taguchi quality loss function. Secondly, taking the overall assembly success rate and the quality loss of the closed loop of each single assembly dimension chain as multiple optimization objectives, improved NSGA-III (NSGA-III-I) is used to solve the simulation data to obtain the Pareto optimal solutions. Finally, considering the production mode of the human-machine intelligent system interaction, fuzzy theory, entropy theory, and the VIKOR method are used to obtain the optimal compromise scheme for selective assembly.
The rest of the paper is organized as follows. In Section 2, the mathematical model and optimization objectives of selective assembly problem for CMPs are introduced. Section 3 presents a case of selective assembly and the assembly of pistons, piston rings, and cylinders. The proposed NSGA-III-I is described in Section 4. In particular, we introduced specific coding, crossover, and mutation strategies to solve the problem of selective assembly. In Section 5, the decision-making method combining fuzzy theory, entropy theory, and VIKOR is described. In Section 6, the NSGA-III-I algorithm is used to optimize the case selective assembly problem and achieved better results compared with other algorithms, and the optimal scheme was obtained successfully from the set of non-dominated solutions by using the proposed decision-making method. Finally, in the last section, conclusions are drawn.

Mathematical Model for Selective Assembly of CMPs
General CMP assembly has the following characteristics: Different quality characteristics (QCs) of the same component constitute different dimension chain closed loops (DCCLs), the same quality characteristics of the same component constitutes a plurality of different DCCLs, and a CMP is composed of many different types and quantities of components [28]. Based on the relationship between DCCLs and QCs, correlation analysis and constraint conditions were established. The optimization mathematical model of the selective assembly of CMPs was established through the assembly success rate and the Taguchi quality loss function.

Correlation Analysis of Components to Be Assembled
In the process of mechanical product assembly, the dimension chain in the assembly link refers to DCCLs composed of QCs from different components in a specific order. DCCLs were used to indicate assembly accuracy during assembly. Therefore, the value of DCCLs represents the actual value of quality requirements (QRs) of the assembly. The correlation between QCs and QRs is the mathematical representation of the assembly structure, and the assembly structure is the basis of selective assembly. The assembly structure relationship of CMPs is shown in Figure 1.
optimize the case selective assembly problem and achieved better results compared with other algorithms, and the optimal scheme was obtained successfully from the set of non-dominated solutions by using the proposed decision-making method. Finally, in the last section, conclusions are drawn.

Mathematical Model for Selective Assembly of CMPs
General CMP assembly has the following characteristics: Different quality characteristics (QCs) of the same component constitute different dimension chain closed loops (DCCLs), the same quality characteristics of the same component constitutes a plurality of different DCCLs, and a CMP is composed of many different types and quantities of components [28]. Based on the relationship between DCCLs and QCs, correlation analysis and constraint conditions were established. The optimization mathematical model of the selective assembly of CMPs was established through the assembly success rate and the Taguchi quality loss function.

Correlation Analysis of Components to Be Assembled
In the process of mechanical product assembly, the dimension chain in the assembly link refers to DCCLs composed of QCs from different components in a specific order. DCCLs were used to indicate assembly accuracy during assembly. Therefore, the value of DCCLs represents the actual value of quality requirements (QRs) of the assembly. The correlation between QCs and QRs is the mathematical representation of the assembly structure, and the assembly structure is the basis of selective assembly. The assembly structure relationship of CMPs is shown in Figure 1. The assembly structure of CMPs is analyzed according to the following steps: Step 1: CMPs are split into individual component, and QC of each component is recorded.
Step 2: Various assembly relations of CMP are found and correspond to QCs of components.
Step 3: Formally express the assembly structure by polychromatic set theory and establish the assembly structure drawing. In this work, multiple dashed block diagrams are used to replace different colors. The assembly structure of CMPs is analyzed according to the following steps: Step 1: CMPs are split into individual component, and QC of each component is recorded.
Step 2: Various assembly relations of CMP are found and correspond to QCs of components.
Step 3: Formally express the assembly structure by polychromatic set theory and establish the assembly structure drawing. In this work, multiple dashed block diagrams are used to replace different colors.

Constraints in Assembly Process
Assembly quality and function of CMPs are determined by the value of each QR. When any QR-value exceeds the upper limit of the objective interval or is lower than the lower limit of the objective interval, the product assembly fails. If all QR-values are within the objective interval, the assembly is said to be qualified. Accordingly, the constraint matrix of QRs is constructed: where y min,n and y max,n denote the lower and upper limits of the design interval value of the nth QR, respectively. Each element in the vector y is the actual value of each QR of the assembly after assembly. When each element in the vector y is within the scope of the upper and lower limits, the assembly is successful. The vector y is as follows:

Taguchi Quality Loss Function
In the process of processing, due to the existence of various random factors, the processed components cannot be completely consistent, and the actual value of the component size has a certain deviation compared with the ideal value. When the deviation is within the acceptable interval, it is called the dimensional tolerance of the component, and the dimensional tolerance of the component determines the QR value of the assembly.
The traditional basis for judging whether a product is qualified is to judge whether the actual QR-value meets the design interval value. As long as it meets the design interval value, the product is considered qualified. Only products that exceed the design interval value are reworked, repaired, or scrapped, and quality loss occurs; according to Taguchi's loss function approach [29], even if it is a qualified product (the actual QR-value is within the design value interval), the fluctuation of the QR-value can still cause losses to users and society. The farther the QR-value is from the ideal value, the greater the loss will be. Therefore, the QR-value should be as close as possible to their ideal value. Two views of quality loss are shown in Figure 2.
lower limit of the objective interval, the product assembly fails. If all QR-values are within the objective interval, the assembly is said to be qualified. Accordingly, the constraint matrix of QRs is constructed: where y min,n and y max,n denote the lower and upper limits of the design interval value of the nth QR, respectively.
Each element in the vector y is the actual value of each QR of the assembly after assembly. When each element in the vector y is within the scope of the upper and lower limits, the assembly is successful. The vector y is as follows:

Taguchi Quality Loss Function
In the process of processing, due to the existence of various random factors, the processed components cannot be completely consistent, and the actual value of the component size has a certain deviation compared with the ideal value. When the deviation is within the acceptable interval, it is called the dimensional tolerance of the component, and the dimensional tolerance of the component determines the QR value of the assembly.
The traditional basis for judging whether a product is qualified is to judge whether the actual QR-value meets the design interval value. As long as it meets the design interval value, the product is considered qualified. Only products that exceed the design interval value are reworked, repaired, or scrapped, and quality loss occurs; according to Taguchi's loss function approach [29], even if it is a qualified product (the actual QR-value is within the design value interval), the fluctuation of the QR-value can still cause losses to users and society. The farther the QR-value is from the ideal value, the greater the loss will be. Therefore, the QR-value should be as close as possible to their ideal value. Two views of quality loss are shown in Figure 2. Taguchi proposed three types of quality loss functions, loss functions for the larger-the-better, the smaller-the-better, the nominal-the-best characteristics, as follows: Taguchi proposed three types of quality loss functions, loss functions for the largerthe-better, the smaller-the-better, the nominal-the-best characteristics, as follows: Larger-the-Better : L(y) = k 1 y 2 Nominal-the-Best : Smaller-the-Better : L(y) = ky 2 (5) where L denotes the Taguchi quality loss, k denotes a constant of proportionality independent of y, y is the actual value of QR, and M is the ideal value of QR. In engineering practice, the nominal-the-best Taguchi quality loss model is often used for the quality loss of each QR of assembly, and it can be guaranteed that the QR-value falls within the design interval [18,30]. For the assembly meeting the design interval, when the QR-value is close to the ideal value, the quality loss caused by the deviation is 0. The quality loss will increase accordingly when the QR-value deviates from the ideal value. If the QR-value exceeds the objective interval value, the quality loss is equal to the cost of product treatment or rework. Therefore, the quality loss of QR is calculated by Formulas (6)-(9): where ∆ denotes the maximum difference between the tolerance boundary and the ideal value, T denotes the width of the design interval, λ denotes the quality loss of assembly failure, and α and β denote the partitioning ratios of the values of the design interval, α + β = 1. The functional image of the nominal-the-best Taguchi quality loss model in assembly is shown in Figure 3.
QR-value falls within the design interval [18,30]. For the assembly meeting th interval, when the QR-value is close to the ideal value, the quality loss caused b viation is 0. The quality loss will increase accordingly when the QR-value devia the ideal value. If the QR-value exceeds the objective interval value, the quali equal to the cost of product treatment or rework. Therefore, the quality loss calculated by Formulas (6)- (9): where ∆ denotes the maximum difference between the tolerance boundary and value, T denotes the width of the design interval, λ denotes the quality loss of a failure, and α and β denote the partitioning ratios of the values of the design inte β = 1. The functional image of the nominal-the-best Taguchi quality loss model i bly is shown in Figure 3.

Many-Objective Optimization Model of Selective Assembly for CMPs
In the selective assembly process, a comprehensive and reasonable evaluati selective assembly scheme of CMPs requires comprehensive consideration of al fiable influencing factors. In selective assembly, the overall success rate and th of each QR close to the ideal value are the most representative evaluation indic scheme.

Many-Objective Optimization Model of Selective Assembly for CMPs
In the selective assembly process, a comprehensive and reasonable evaluation of the selective assembly scheme of CMPs requires comprehensive consideration of all quantifiable influencing factors. In selective assembly, the overall success rate and the degree of each QR close to the ideal value are the most representative evaluation indices of the scheme.

Assessment Based on Assembly Success Rate
In the selective assembly of existing components to be assembled, more successful assemblies should be obtained as far as possible under the condition of meeting assembly requirements, and the surplus of mismatched components should be reduced. Assembly success rate is defined in the range [0, 1]. The assembly success rate in optional scheme X can be calculated by Formula (10): where n a denotes the number of assemblies to be assembled; n s denotes number of qualified assemblies obtained through selective assembly; X denotes the scheme for selective assembly; and η s denotes the success rate of selective assembly.

Assessment Based on Quality Loss of Individual QR
In previous studies, the multidimensional chain was often converted into a singledimensional chain when the quality loss of the assembly was modeled, and the simple linear addition of the quality loss was used to set an optimization objective. This will lead to product quality requirement confusion, and it is not easy to judge the effect of the assembly scheme. The optimization effect of product assembly is not ideal, and the performance of schemes with the same quality loss varies greatly. Based on the characteristics of general CMPs, the proposed work takes the average quality loss function of a single QR as the optimization objective function as calculated by Formulas (11) and (12): where T i denotes the width of the design interval for the ith QR, y ij denotes the actual value of the ith QR of the jth assembly, and M i denotes the objective value of the ith QR. L(y ij ) denotes the quality loss of the ith QR of the jth assembly, and L(y i ) denotes the mean value of the sum of the quality loss of the ith QR of all assemblies in the selective assembly scheme; q denotes the number of components of the same type.

Many-Objective Optimization Model of Selective Assembly for CMPs
The pros and cons of CMP batch selective assembly schemes are jointly evaluated by the overall assembly success rate and product quality loss. To unify the optimization direction of a set of objective functions to the minimum direction, the first objective function is deformed. From the formula in this chapter, the mathematical model (Formula (13)) of the many-objective optimization for selective assembly can be obtained as follows: where m denotes the number of optimization objective functions, m = n + 1, n denotes the number of QR; f (X) denotes the objective vector, f 1 (X), f 2 (X), . . . , f m (X) denote each element value in the vector, and X denotes the selective assembly scheme; x p denotes the full permutation of the pth type of components; O q denotes the domain of x 1 , x 2 , . . . , x p , which is the full permutation of 1 to q; q denote the number of components of the same type.

Case Analysis
An example of the selective assembly of pistons and cylinders proposed by Asha [20] is used to verify the selective assembly of CMPs in this work. The assembly structure is shown in Figure 4.

Case Analysis
An example of the selective assembly of pistons and cylinders proposed [20] is used to verify the selective assembly of CMPs in this work. The assemb ture is shown in Figure 4.

Components QCs and QRs of the Assembly
The assembly is composed of three components: a piston, a piston ring, an inder. The considered QCs of the components are listed in Table 1.

Components QCs and QRs of the Assembly
The assembly is composed of three components: a piston, a piston ring, and a cylinder. The considered QCs of the components are listed in Table 1. The QCs of all the components of the piston and cylinder assembly are analyzed and QRs that need to be analyzed and calculated during the selective assembly process of the assembly are obtained, as shown in Figure 5. The QCs of all the components of the piston and cylinder assembly are analyzed and QRs that need to be analyzed and calculated during the selective assembly process of the assembly are obtained, as shown in Figure 5.    If the components are assembled through the principle of interchangeability, the tolerances of four QRs will be the sum of the tolerances of the corresponding component QCs, which is the maximum tolerance. Consider the center value of each design tolerance interval as the ideal value of each QR, and the design tolerance intervals of each QR of the assembly are listed in Table 2. According to the design tolerance interval of each QR of the assembly, the constraint matrix is obtained as follows: The ideal value vector of QRs is as follows: The objective tolerance interval of QRs is as follows:

Optimization Objectives of the Case Selective Assembly
According to the many-objective selective assembly optimization model established in this work and the assembly QRs of the piston and cylinder assembly in the case study, the optimization objectives in the selective assembly process of the piston and cylinder group are established as follows:

Data Simulation of Each Component
In this work, the purpose of the proposed method was to maximize the success rate of assembly while reducing the quality loss of each quality requirement as much as possible. To verify the feasibility of the scheme, the Monte Carlo simulation method was used to obtain the QC data of each component. A total of 300 sets of components were simulated and randomly divided into 6 groups, each with 50 components. Simulating the component QC data in the actual production process through the discrete normal distribution probability mass function (Formula (14)) is defined as: where H denotes the normal distribution of QCs of the component, Z denotes the set of integers, and the size table of the components is listed in Table A1.

NSGA-III-I
The selective assembly model of the CMPs constructed by the above content was a many-objective optimization model with multiple objectives that needed to be optimized at the same time, and the optimization objectives were often greater than three. An NSGA-III-I selective assembly many-objective optimization algorithm suitable for this model was designed.
The NSGA-III is a new many-objective optimization algorithm proposed on the basis of the NSGA-II [31]. The NSGA-II selects individuals of the same non-dominant level by calculating the distance of individual crowding degrees. It was believed that the larger the crowding degree distance is, the better. This method was particularly suitable for solving the optimization problem of two objectives. However, when it solves optimization problems with more than three objectives, the individuals selected by crowding distance sorting are not uniformly distributed on the non-dominated layer, which easily falls into the local optimum, which affects the performance of the algorithm. The NSGA-III selects individuals based on the method of reference points and optimizes the selection by calculating the shortest distance between individuals and reference points. When facing many-objective optimization problems with three or more objectives, the algorithm can be ensured to converge to a uniformly distributed Pareto surface [32].
The NSGA-III algorithm is limited by the coding method and completely random initial solution generation rules, resulting in the algorithm performance being greatly dependent on the initial solution set, and the feasible solution space of the constructed selective assembly model increases exponentially with the increase in components. It greatly increases the risk of obtaining a local optimal solution. To improve the global searching ability of the NSGA-III, a new mapping method between genotype and phenotype was selected, and the idea of a simulated annealing algorithm was introduced into crossover and mutation operators to enhance the ability of the algorithm to obtain the global Pareto front (PF).

Pareto Optimal Solution
The solution set finally obtained by the NSGA-III algorithm is the Pareto optimal solution set. When the optimization direction is the minimum value direction, the solution set is defined as follows: For any two decision variables X a , X b ∈ X f , where X f is the feasible solution set: If and only if the above conditions are met, X b is dominated by X a . In the solution space, there is a solution set X, which is not dominated by any solution; it is the Pareto optimal solution [33].

Input Module
For the experimental problems considered in this paper, the required input data are as follows: • Problem data: (1) Batch size; (2) Number of component types in the assembly; Dimensional data of components.

Encoding
Selective assembly is a many-objective combinatorial optimization problem. The way to solve the problem is to find the best possible selective assembly scheme. Traditional binary coding will produce a large number of infeasible solutions when contending with selective assembly problems, which increases the difficulty of the algorithm in solving the problem. Kannan et al. [5] proposed a classic coding method to solve the selective assembly problem using genetic algorithms ( Figure 6a). However, this kind of coding scheme still encounters failure in achieving a one-to-one correspondence between the selective assembly scheme and the chromosome. In the case of different genotypes, it can match the same selective assembly scheme, which will cause the performance of the algorithm to decrease. Therefore, Rezaei Aderiani et al. [25] proposed a new coding scheme that can ensure that the dyeing and selection scheme is one-to-one, and the dotted lines indicate that the genes are fixed and immutable, as shown in Figure 6b. Taking an assembly consisting of two components A and B as an example, five for each component, (A1B4), (A2B3), (A3B5), (A4B1), and (A5B2), this combination can be expressed as:

Encoding
Selective assembly is a many-objective combinatorial optimization problem. The way to solve the problem is to find the best possible selective assembly scheme. Traditional binary coding will produce a large number of infeasible solutions when contending with selective assembly problems, which increases the difficulty of the algorithm in solving the problem. Kannan et al. [5] proposed a classic coding method to solve the selective assembly problem using genetic algorithms ( Figure 6a). However, this kind of coding scheme still encounters failure in achieving a one-to-one correspondence between the selective assembly scheme and the chromosome. In the case of different genotypes, it can match the same selective assembly scheme, which will cause the performance of the algorithm to decrease. Therefore, Rezaei Aderiani et al. [25] proposed a new coding scheme that can ensure that the dyeing and selection scheme is one-to-one, and the dotted lines indicate that the genes are fixed and immutable, as shown in Figure 6b. Taking an assembly consisting of two components A and B as an example, five for each component, (A1B4), (A2B3), (A3B5), (A4B1), and (A5B2), this combination can be expressed as:

Initialization of Algorithm Parameters
Initialize the population and randomly generate N chromosomes; Initialize the number of iterations, Gen = 0; Initialize to reference points.

Initialization of Algorithm Parameters
Initialize the population and randomly generate N chromosomes; Initialize the number of iterations, Gen = 0; Initialize to reference points.

New Population Generation Module
The new population generation module was divided into a crossover module and a mutation module. In NSGA-III, the gene mutation method with fixed crossover and mutation probability was adopted. This strategy keeps the crossover and mutation probability unchanged during the operation of the algorithm and cannot consider the global search performance and convergence performance required by the algorithm. Aiming at the many-objective optimization problem in the CMP selective assembly, the adaptive change idea of a simulated annealing algorithm was used to improve crossover and mutation operators to enhance the searching ability and convergence of the algorithm.

New Crossover Operator
Notably, when the crossover probability of the algorithm was small, it was beneficial to maintaining population diversity and increasing global search ability. In contrast, when the crossover probability was large, the efficiency of local search was strong, and the convergence ability of the algorithm was improved. Therefore, the crossover operator was improved according to the crossover probability characteristics of the algorithm. In particular, the crossover probability was small in the early stage of the algorithm operation to maintain population diversity, and the crossover probability was large in the late stage to improve the convergence ability of the algorithm. The crossover strategy (Formula (15)) is as follows: where ProC denotes the crossover probability, p C1 and p C2 denote the crossover parameters, and Gen denotes the number of iterations. In this work, p C1 = 0.6 and p C2 = 50 are selected.
Step 1: Select crossed chromosomes and generate a random number r for each pair of chromosomes, if r < ProC, corresponding chromosomes are selected for the crossover operation.
Step 2: The chromosomes selected for crossover are paired, and crossover operations are performed between the same substrings of different chromosomes. To avoid invalid chromosomes, the crossover point is placed between the second and thirds substrings, and a single point crossover is performed. An example of a two-chromosome crossing operation to produce two offspring chromosomes is shown in Figure 7.
Processes 2021, 9, x FOR PEER REVIEW 13 of 34 particular, the crossover probability was small in the early stage of the algorithm operation to maintain population diversity, and the crossover probability was large in the late stage to improve the convergence ability of the algorithm. The crossover strategy (Formula (15)) is as follows: where ProC denotes the crossover probability, p C1 and p C2 denote the crossover parameters, and Gen denotes the number of iterations. In this work, p C1 = 0.6 and p C2 = 50 are selected.
Step 1: Select crossed chromosomes and generate a random number r for each pair of chromosomes, if r < ProC, corresponding chromosomes are selected for the crossover operation.
Step 2: The chromosomes selected for crossover are paired, and crossover operations are performed between the same substrings of different chromosomes. To avoid invalid chromosomes, the crossover point is placed between the second and thirds substrings, and a single point crossover is performed. An example of a two-chromosome crossing operation to produce two offspring chromosomes is shown in Figure 7.

New Mutation Operator
The mutation operator locally changes the genes of the chromosomes, hoping to create better chromosomes. When the algorithm mutation probability is too low, it will greatly increase the risk of falling into the local optimum; when the algorithm mutation probability is too high, it will cause the algorithm solution set to fail to converge. Therefore, by improving the mutation strategy of the algorithm, this study makes the mutation probability relatively high at the beginning of the algorithm, increases the diversity and richness of the population, and gradually decreases with the increase in the number of

New Mutation Operator
The mutation operator locally changes the genes of the chromosomes, hoping to create better chromosomes. When the algorithm mutation probability is too low, it will greatly increase the risk of falling into the local optimum; when the algorithm mutation probability is too high, it will cause the algorithm solution set to fail to converge. Therefore, by improving the mutation strategy of the algorithm, this study makes the mutation probability relatively high at the beginning of the algorithm, increases the diversity and richness of the population, and gradually decreases with the increase in the number of iterations, making it easier to retain more excellent individuals and improving the convergence performance of the algorithm. The mutation strategy (Formula (16)) is as follows: where ProM denotes the mutation probability, p m1 and p m2 denote mutation parameters, and Gen denotes the number of iterations. In this work, p m1 = 0.6 and p m2 = 50 are selected. Mutation operations are used to target genes on each chromosome. To avoid the occurrence of invalid chromosomes, the neighborhood swapping method is used to generate a random number r for each gene; if r < ProM during the mutation process, the selected gene is swapped with the previous gene. An example of the mutation operation is shown in Figure 8. selected gene is swapped with the previous gene. An example of the mutation operation is shown in Figure 8.

Merging of Parent and Offspring Populations
The combined population R t with population size 2N was obtained by merging the parent population P t with population size N and the offspring population Q t with population size N. R t is the transition population between population P t and the next population P t+1 .

Evaluation Module and Environment Selection Module
The fitness vector elements of each chromosome correspond to the five objective function values of f 1 (X), f 2 (X), f 3 (X), f 4 (X), and f 5 (X). The five objective function values of N chromosomes are calculated.

Selected Individuals
The transition group R t with a population size of 2N is merged by the current group P t and its progeny group Q t . To make the size of P t+1 meet the requirements of population size and ensure that the overall performance is better than P t , a non-dominated sorting operation is performed on R t , and individuals are selected with higher non-dominated levels and added to P t+1 until the size of P t+1 cannot accommodate individuals at the next level of dominance, and this dominance level is recorded as F l .
The associated reference point operations on individuals in the current population (excluding individuals in F l ) are performed, the reference point j with the least associated individuals is selected and its associated number is recorded as ρ j . Then, the situation of this reference point can be elaborated as follows:

Merging of Parent and Offspring Populations
The combined population R t with population size 2N was obtained by merging the parent population P t with population size N and the offspring population Q t with population size N. R t is the transition population between population P t and the next population P t+1 .

Evaluation Module and Environment Selection Module
The fitness vector elements of each chromosome correspond to the five objective function values of f 1 (X), f 2 (X), f 3 (X), f 4 (X), and f 5 (X). The five objective function values of N chromosomes are calculated.

Selected Individuals
The transition group R t with a population size of 2N is merged by the current group P t and its progeny group Q t . To make the size of P t+1 meet the requirements of population size and ensure that the overall performance is better than P t , a non-dominated sorting operation is performed on R t , and individuals are selected with higher non-dominated levels and added to P t+1 until the size of P t+1 cannot accommodate individuals at the next level of dominance, and this dominance level is recorded as F l .
The associated reference point operations on individuals in the current population (excluding individuals in F l ) are performed, the reference point j with the least associated individuals is selected and its associated number is recorded as ρ j . Then, the situation of this reference point can be elaborated as follows: (1) If there is no individual associated with the reference point in F l , a new reference point vector is selected. (2) If the number of population individuals associated with this reference point is 0 (ρ j = 0), but there are individuals in F l that are related to this reference point vector, then the individual with the smallest distance from it is found and extracted from F l , which is then added it to the selected next-generation population, set ρ j = ρ j + 1. (3) If ρ j > 0, there are multiple individuals in F l that are associated with this reference point vector, and then individuals in F l are randomly selected that are associated with the reference point until the population size is N.

Associated Reference Point
According to the fitness vector of all individuals in the population, the minimum value of each objective function is calculated, the ideal point is updated, the hyperplane is updated, each individual is associated to the reference point and reference vector according to the strategy, and the number of individuals associated with the reference point is calculated. The specific association strategy is as follows: Step 1 Generate hyperplane: First, select the minimum value of each objective function in the current population, calculate the ideal point, take the ideal point as the origin and the objective function as the coordinate axis, and normalize the individual using Formula (17): where z min j denotes the ideal point; f j (X) denotes the jth fitness value of each individual in the population; and f j (X) denotes jth fitness value of each individual after normalization in the population.
Next, the achievement scalarizing function (ASF(X , w) = m max j = 1 f j (X)/ w j ) is used to calculate the extreme points. An extreme point refers to the individual whose objective function value of one dimension is very large and the objective function value of the other dimensions is very small; then, the population is traversed, and the extreme points z max i of each dimension are found. Then, a hyperplane is constructed based on these extreme points, and the intersection of the hyperplane and the coordinate axis is the intercept a i . Taking the three-object problem as an example, the hyperplane is obtained as shown in Figure 9. is updated, each individual is associated to the reference point and reference vector according to the strategy, and the number of individuals associated with the reference point is calculated. The specific association strategy is as follows: Step 1 Generate hyperplane: First, select the minimum value of each objective function in the current population, calculate the ideal point, take the ideal point as the origin and the objective function as the coordinate axis, and normalize the individual using Formula (17): where z j min denotes the ideal point; f j (X) denotes the jth fitness value of each individual in the population; and f j ' (X) denotes jth fitness value of each individual after normalization in the population.
Next, the achievement scalarizing function (ASF(X , w) = m max j=1 f j ' (X) / w j ) is used to calculate the extreme points. An extreme point refers to the individual whose objective function value of one dimension is very large and the objective function value of the other dimensions is very small; then, the population is traversed, and the extreme points z i max of each dimension are found. Then, a hyperplane is constructed based on these extreme points, and the intersection of the hyperplane and the coordinate axis is the intercept a i . Taking the three-object problem as an example, the hyperplane is obtained as shown in Figure 9. In NSGA-III, it is necessary to obtain a hyperplane equidistant from the origin and normalize the individuals of the population according to the intercept of the coordinate axes so that a standard hyperplane can be obtained in each generation. After this transformation operation, the individuals of the population are normalized to the plane where the reference point is. The conversion Formula (18) is as follows: In NSGA-III, it is necessary to obtain a hyperplane equidistant from the origin and normalize the individuals of the population according to the intercept of the coordinate axes so that a standard hyperplane can be obtained in each generation. After this transformation operation, the individuals of the population are normalized to the plane where the reference point is. The conversion Formula (18) is as follows: Step 2 Generate reference points: This study adopted the reference point establishment method proposed by Das and Dennis [34]: The reference point is placed on a normalized m − 1 dimensional standard hyperplane, which is evenly inclined to all objective axes and has an intercept on each axis. H points are uniformly generated on the hyperplane, and the boundary crossing method is used to construct weights. If the objective dimension is divided into p equal parts, the total number of reference points H in the m-objective problem is ( m + p − 1 p ).
Taking the three-objective problem as an example, the reference point is created in the triangle plane with (0, 0, 1), (0, 1, 0), and (1, 0, 0) as the vertices. If the coordinate axis is divided into four equal parts, p = 4, H = ( 3 + 4 − 1 4 )= 15. The distribution of reference points is shown in Figure 10. This study adopted the reference point establishment method proposed by Das and Dennis [34]: The reference point is placed on a normalized m−1 dimensional standard hyperplane, which is evenly inclined to all objective axes and has an intercept on each axis. H points are uniformly generated on the hyperplane, and the boundary crossing method is used to construct weights. If the objective dimension is divided into p equal parts, the total number of reference points H in the m-objective problem is m + p 1 p .
Step 3 Associated reference point: Each reference vector is constructed by connecting each reference point and the origin. For each individual in the population, the reference vector was traversed to find the reference point nearest to the individual, and the association information and distance were recorded. The association between the individual and the reference point was completed. The distance between the individual and the reference vector will be described by the Euclidean distance.

Termination Criterion
When the number of iterations reaches the maximum number or other relevant setting conditions are met, the algorithm terminates. In the experiment in this paper, the maximum number of iterations of the algorithm is set to 500.

Output Module
The output of the algorithm is the 500th generation population, and the population size is N. Decoding can obtain N non-dominated selective assembly schemes.
Step 3 Associated reference point: Each reference vector is constructed by connecting each reference point and the origin. For each individual in the population, the reference vector was traversed to find the reference point nearest to the individual, and the association information and distance were recorded. The association between the individual and the reference point was completed. The distance between the individual and the reference vector will be described by the Euclidean distance.

Termination Criterion
When the number of iterations reaches the maximum number or other relevant setting conditions are met, the algorithm terminates. In the experiment in this paper, the maximum number of iterations of the algorithm is set to 500.

Output Module
The output of the algorithm is the 500th generation population, and the population size is N. Decoding can obtain N non-dominated selective assembly schemes.
The basic search process of the proposed algorithm is shown in Figure 11.
Processes 2021, 9, x FOR PEER REVIEW 17 of 34 Figure 11. Flowchart of the proposed NSGA-III-I for batch selective assembly.

VIKOR Method to Sort the Solution Set
The Pareto optimal solution set obtained by the NSGA-III-I algorithm solving the selective assembly high-dimensional many-objective optimization model is actually a non-dominated solution set, which includes a large number of mutually non-dominated solutions. It is impossible to determine the priority of each scheme from the numerical relationship, so the schemes of selective assembly need to be sorted and selected by means of human-machine intelligent system interaction. The VIKOR method is a multi-criteria decision-making method proposed by Opricovic [35], which is an optimal Figure 11. Flowchart of the proposed NSGA-III-I for batch selective assembly.

VIKOR Method to Sort the Solution Set
The Pareto optimal solution set obtained by the NSGA-III-I algorithm solving the selective assembly high-dimensional many-objective optimization model is actually a non-dominated solution set, which includes a large number of mutually non-dominated solutions. It is impossible to determine the priority of each scheme from the numerical relationship, so the schemes of selective assembly need to be sorted and selected by means of human-machine intelligent system interaction. The VIKOR method is a multi-criteria decision-making method proposed by Opricovic [35], which is an optimal compromise decision-making method based on ideal points, as well as an optimal compromise solution. It is developed from the L p-metric distance measure function (Formula (19)): where X i denotes the ith scheme; p denotes the distance parameter of the L p-metric distance measure function; f ij denotes the evaluation value of the ith scheme in the jth index; f * j and f − j denote the best and worst evaluation values of the jth index; ω j denotes the weight of the jth index.

Calculation of the Weight of Each Indicator
In the decision-making process of human-machine intelligent system interaction, it is necessary to assign weight to each optimization objective, and the index weight is a measure reflecting the relative importance of each optimization objective. In this paper, the entropy weight method is used to calculate the objective weight of the index. The triangular fuzzy number is used to fuzzily quantify the engineering decision maker language variables, and then the subjective weight of the indicators is solved. Finally, the subjective weight and the objective weight are combined into a combination weight with high robustness.

Entropy Method to Determine the Objective Weight
The entropy method is an approach to obtain the weight [36]. According to the basic principles of information theory, information is a measure of the order degree of a system, and entropy is a measure of the disorder degree of a system. If the information entropy of the indicator is smaller, then the amount of information provided by the indicator is larger, and it should play a larger role in the comprehensive evaluation such that the weight should be higher.
According to the definition of information entropy in information theory, the information entropy of a set of data, and the objective weights of each indicator are calculated using Formulas (20)- (22): where E j denotes the information entropy, p ij denotes the proportion of f ij in the jth indicator data, n is the number of decision-making schemes, and f ij denotes the value of the jth objective function in the ith scheme after normalization. ω j denotes the objective weight of the jth indicator.

Triangular Fuzzy Numbers to Determine Subjective Weights
Taking into account the uncertainty and ambiguity in the decision making process, a triangular fuzzy number was introduced to establish a subjective weight matrix [37]. The correspondence between linguistic variables and triangular fuzzy numbers is determined, as listed in Table 3: Different engineering decision makers may have different views on each index. The fuzzy weight vector of decision makers B k on the optimization objective is ω jk = (ω jk1 , ω jk2 , ω jk3 )| j = 1, 2, . . . , m . Therefore, after gathering different decision makers to evaluate the indicators, each element of the subjective weight matrix is expressed using Equations (23) and (24): Then, the subjective weight of the indicator is obtained. Because each element of the subjective weight matrix is a triangular fuzzy number, the mean area method is used for defuzzification to obtain the subjective weight of each indicator (Equation (25)):

Determination of the Combination Weight
In order to balance the impact of the subjective and objective factors on the result of the selective assembly scheme, multiplication and division were selected to determine the combined weight of the subjective weight and the objective weight. The calculation method of the combined weight uses Equation (26): where γ and γ denote the proportion coefficients of the objective weight and subjective weight in the total weight, respectively. Generally, γ +γ = 1. In this work, γ = γ = 0.5.

VIKOR Method for Multicriteria Decision-Making
VIKOR's basic point of view is to fully consider the subjective preferences of decision makers and, at the same time, to trade-off the limited decision-making schemes by maximizing group benefits and minimizing individual regrets, making the decision results more reliable. The compromise solution F c is the closest to the ideal solution F * , which is the result of mutual concessions between the two criteria, as shown in Figure 12. (3) Calculate the group utility value S i and individual regret value R i of each Pareto optimal solution by using Equations (27) and (28): (4) Calculate the compromise value Q i of each Pareto optimal solution scheme by using Equation (29): where S Max , S Min , R Max , and R Min are the maximum and minimum values of S i and R i , respectively. (5) Sort the Pareto optimal solution set X according to the increasing Q i -value: X 1 , X 2 , …, X i , X n . If X 1 is the solution with the smallest Q-value and both Condition 1 and Condition 2 are satisfied, X 1 is the stable optimal solution. Condition 1: Q(X 2 ) − Q(X 1 ) ≥1/(n − 1); Condition 2: After the solutions are sorted based on the Q-value, the S-value, or R-value of the first-ranked solution X 1 is better than the second-ranked solution X 2 .
If the above two conditions cannot be established at the same time, a compromise solution will be obtained, which can be divided into two cases: Case 1: If only Condition 2 is not satisfied, the compromise solution is X 1 , X 2 ; Case 2: As long as Condition 1 is not satisfied, the compromise solution is X 1 , X 2 , …, X i , where i is the maximized value determined by Q(X i ) − Q(X 1 ) ≥1/(n − 1).  (27) and (28): (4) Calculate the compromise value Q i of each Pareto optimal solution scheme by using Equation (29): where S max , S min , R max , and R min are the maximum and minimum values of S i and R i , respectively. (5) Sort the Pareto optimal solution set X according to the increasing Q i -value: X 1 , X 2 , . . ., X i , X n . If X 1 is the solution with the smallest Q-value and both Condition 1 and Condition 2 are satisfied, X 1 is the stable optimal solution.
Condition 1: Q(X 2 ) − Q(X 1 ) ≥ 1/(n − 1); Condition 2: After the solutions are sorted based on the Q-value, the S-value, or R-value of the first-ranked solution X 1 is better than the second-ranked solution X 2 .
If the above two conditions cannot be established at the same time, a compromise solution will be obtained, which can be divided into two cases: Case 1: If only Condition 2 is not satisfied, the compromise solution is X 1 , X 2 ; Case 2: As long as Condition 1 is not satisfied, the compromise solution is X 1 , X 2 , . . ., X i , where i is the maximized value determined by Q(X i ) − Q(X 1 ) ≥ 1/(n − 1).

Analysis of Simulation Results of Different Algorithms
In order to verify the effectiveness of the NSGA-III-I algorithm in solving manyobjective selective assembly optimization problems, a comparative simulation experiment was designed for the NSGA-III-I algorithm, interchangeable assembly (IA), NSGA-II, and NSGA-III. The basic parameters and operating environment of each algorithm were consistent. The original code for each algorithm comes from PlatEMO [38]. All of the algorithms were programmed with MATLAB and conducted on an AMD Ryzen 5 2400 G Processor computer. The parameters of each algorithm are as shown in Table 4. The hypervolume (HV) index was used to compare the performance of different algorithms for selective assembly problems. Suppose the solution set obtained by the algorithm is X p , Ref = (r 1 , r 2 , · · · , r m ) is the reference point and is dominated by any of the solutions in X p . The HV value is the spatial hypervolume of the solution set X p bounded by the reference point Ref. It is calculated by using Formula (30): where Leb(X p ) denotes the Lebesgue measure of the solution set X p and [ f 1 , r 1 ] × [ f 2 , r 2 ] × · · · × [ f m , r m ] denote a hypercube surrounded by all points that are dominated by the solution in X p but not dominated by Ref.
The HV index is sensitive to the selection of reference point Ref when measuring the performance of different algorithms. To minimize the influence of the selection of reference points on the HV value, the final solution set of multiple algorithms was combined, and nondominated sorting was conducted to obtain the approximate real PF for the same set of experimental data, taking the reference point Ref = (1, 1, · · · , 1) [38]. The HV index is a positive index, and the larger the HV value is, the better the comprehensive performance of the algorithm is. When the HV value is 0, it indicates that the comprehensive performance of the algorithm is poor and that the solution set does not converge to the PF.
The 300 sets of component data were equally divided into 6 batches, each with 50 sets of components. Each algorithm was run 6 times per batch to reduce the influence of randomness, and the performance of each algorithm in the selection of assembly problems was compared according to the running results. The optimal value of the HV value of each algorithm solution set, the median value, the worst value and the mean value are listed in Table 5. Rows 1-4 in each group show the optimal, median, worst, and average HV values of the IA, NSGA-II, NSGA-III, and NSGA-III-I, respectively. The bolded data are the maximum value of each row. NSGA-III-I obtained the maximum optimal value, maximum worst value, maximum median value and maximum mean value of the HV value in six groups of simulation experiments. NSGA-III-I has better overall performance compared with other algorithms.  The content in bold is the maximum value of each row.
To reflect the advantages and disadvantages of the Pareto optimal solution set of each algorithm more clearly, the Pareto optimal solution set obtained by the group of experiments of the first batch (batch size = 50) was used to describe the algorithms. Six random simulation experiments of each algorithm used the HV value to sort, and then the solution sets of different algorithms with the same serial number were classified and compared. The parallel coordinate diagrams of the solution sets of the three algorithms in the first group of experiments are shown in Figure 13.
In this group of experiments, a partial dominance relationship existed in the corresponding solution set after sorting the HV values of the solution set of each algorithm. The average value of each objective of the IA method was selected, and the partial nondominated solutions obtained by each algorithm were compared, as shown in Table A2.
In the six repeated experiments of the first batch, under the condition that the solution sets obtained by each algorithm are the current non-dominated solutions, the solution sets obtained by NSGA-III-I in the six experiments all have a certain number of solutions dominating the solutions of other methods. It shows that the NSGA-III-I has better performance in solving the problem of many-objective optimization in the selective assembly than other methods in the study, and the success rate increased from 57% (IA method) to 96-98% (NSGA-III-I) while the quality loss of each quality requirement was significantly reduced. In this group of experiments, a partial dominance relationship existed in the c sponding solution set after sorting the HV values of the solution set of each algor The average value of each objective of the IA method was selected, and the p non-dominated solutions obtained by each algorithm were compared, as shown in A2.
In the six repeated experiments of the first batch, under the condition that the tion sets obtained by each algorithm are the current non-dominated solutions, the tion sets obtained by NSGA-III-I in the six experiments all have a certain number lutions dominating the solutions of other methods. It shows that the NSGA-III-I ha ter performance in solving the problem of many-objective optimization in the sele assembly than other methods in the study, and the success rate increased from 57% method) to 96-98% (NSGA-III-I) while the quality loss of each quality requiremen significantly reduced.

Analysis of Convergence of Each Algorithm
In the experiments, the convergence process of each algorithm is described b variation of the minimum value of each objective function in each generation. The value of six repeated experiments is used to reduce the error (Figure 14).

Analysis of Convergence of Each Algorithm
In the experiments, the convergence process of each algorithm is described by the variation of the minimum value of each objective function in each generation. The mean value of six repeated experiments is used to reduce the error (Figure 14). The results show that each objective function value of the NSGA-II and NSGA-III decreases rapidly in the iterative process but gradually falls into the local optimum around 200 generations, and it is difficult to jump out. In the NSGA-III-I, each objective function value keeps a relatively stable decline rate and can finally explore a solution set closer to the real Pareto frontier.

Analysis of the Influence of Components Batch Size on Algorithm
To further study the influence of batch size and quantity on the performance of each algorithm, the same basic parameters and operating environment were used to repeat the test on the same 300 sets of components. The experiment is repeated for the same 300 mating components with 12 batches of 25 components each, 4 batches of 75 components each, and 3 batches of 100 components each. The results are shown in Tables A3-A5, re- The results show that each objective function value of the NSGA-II and NSGA-III decreases rapidly in the iterative process but gradually falls into the local optimum around 200 generations, and it is difficult to jump out. In the NSGA-III-I, each objective function value keeps a relatively stable decline rate and can finally explore a solution set closer to the real Pareto frontier.

Analysis of the Influence of Components Batch Size on Algorithm
To further study the influence of batch size and quantity on the performance of each algorithm, the same basic parameters and operating environment were used to repeat the test on the same 300 sets of components. The experiment is repeated for the same 300 mating components with 12 batches of 25 components each, 4 batches of 75 components each, and 3 batches of 100 components each. The results are shown in Tables A3-A5 respectively.
The results show that when the number of components in each batch is small (batch size = 25), the proposed algorithm obtains the maximum value in the optimal value, the maximum value in the median value and the maximum value in the worst value of HV in most experimental groups. With the increase in the number of components in each batch, the algorithm obtained the maximum value of the optimal value, the maximum value of the median value, and the maximum value of the worst value of all experimental groups.
As the number of components in each batch increases, some algorithms have gradually failed to converge ( Figure 15). The results show that the IA method completely fails to converge to the approximate PF. With the increase in the number of components in each batch, both the NSGA-II and NSGA-III fail to converge to varying degrees. Only the NSGA-III-I algorithm can still converge to the approximate PF, with better performance than other algorithms and a better selective assembly scheme obtained.

The VIKOR Method to Sort the Optimal Compromise Scheme
Taking the first set of experiments with 50 components per batch as an example, the Pareto solution set obtained by the NSGA-III-I algorithm is used as a decision matrix. The ideal solution in this experiment F -= Min f(x) = (0.0200, 0.1658, 0.2320, 0.1731, 0.1537). The step-by-step procedures and results are given below.
Step 1: The entropy method to obtain the objective weight (listed in Table 6). Step 2: Triangular fuzzy numbers to determine the subjective weights of the assem-

The VIKOR Method to Sort the Optimal Compromise Scheme
Taking the first set of experiments with 50 components per batch as an example, the Pareto solution set obtained by the NSGA-III-I algorithm is used as a decision matrix. The ideal solution in this experiment F − = Min f (x)= (0.0200, 0.1658, 0.2320, 0.1731, 0.1537). The step-by-step procedures and results are given below.
Step 1: The entropy method to obtain the objective weight (listed in Table 6).
Step 2: Triangular fuzzy numbers to determine the subjective weights of the assembly quality requirements. Six engineering decision makers were invited to evaluate the importance of each optimization objective in the cylinder-piston assembly selection optimization model (listed in Table 7). The comparison table converts semantic variables into triangular fuzzy numbers (listed in Table 8). The defuzzification method (Equation (24)) is performed to obtainthe subjective weights (listed in Table 9).  Step 3: Combination weightings are calculated by using Equation (25) (listed in Table 10). Step 4: Calculation of the S-value, R-value, and Q-value of each Pareto solution in turn. The six optimal schemes sorted according to the Q-values are listed in Table 11. Table 11. Information on the 6 schemes with the best Q-value sorting.

Experiment No
Chromosome No. Step 5: Discussions and decision making based on the above results. The analysis of the calculation results in Table 11 shows that (1) Q(X 2 ) − Q(X 1 ) ≥ 1/(126 − 1), and (2) when sorting according to the Q-value, scheme X 1 is the optimal choice; when sorting according to the R value, scheme X 1 is still the optimal choice; and when sorting according to the S value, scheme X 1 is superior to scheme X 2 . According to VIKOR theory, when the decision maker has a compromising attitude towards all schemes (v = 0.5), scheme X 1 can be selected as the optimal scheme among all non-dominated schemes (listed in Table 12).
According to the above five steps, we conclude that scheme X 1 is the optimal compromise solution. Table 12 shows the detailed content of scheme X 1 , where each row represents the sequence number arrangement of the same type of components. The 50 pistons, 50 piston rings, and 50 cylinders are assembled in one-to-one correspondence. For example, piston #1 lines up with piston ring #18 and cylinder #35, and they are all located at the first position in their respective arrangement sequence. In particular, piston #44, piston ring #19, and cylinder #23 is the only set that cannot be successfully assembled in the 50 sets of components.

Conclusions
This research focuses on the optimization of selective assembly of complex mechanical products (CMPs) under intelligent system environment conditions. A general method based on the NSGA-III-I and VIKOR was proposed to solve the many-objective optimization and decision-making problem in the selective assembly of CMPs, which could optimize the assembly effect more effectively and select the assembly scheme rationally. The main contributions are as follows: (1) A mathematical model of selective assembly based on the success rate and the Taguchi quality loss was constructed, dimensional constraints and quality requirements were established through assembly structure analysis, and multiple optimization objectives s were established on this basis.
(2) For several optimization objectives, the NSGA-III-I was proposed, and the simulation case was solved. Experiments were designed to compare the NSGA-III-I with the IA, NSGA-II, and NSGA-III methods, and the performance changes of each algorithm under different batch conditions were studied. The results show that the proposed algorithm has obvious advantages in solving many-objective selective assembly problems. The proposed method increases the assembly success rate from 57% (IA) to 96-98% (NSGA-III-I), while reducing quality loss significantly. It effectively avoids falling into the local optimum compared with the NSGA-II and NSGA-III. (3) Considering the production mode of human-machine intelligent system interaction, the optimal compromise solution is obtained by using fuzzy theory, entropy theory, and the VIKOR method. The planning scheme can simultaneously maximize group benefits and minimize individual regret.
Future research work will add a quantitative method and constraint conditions of form and position tolerance on the basis of the proposed selective assembly model so that the method can be more widely used in the selective assembly process of almost all components. Therefore, it can be integrated into a digital manufacturing execution system (MES) and systematically applied to CMP manufacturing enterprises. After designing for specific problems, the proposed method can be applied to many different types of many-objective combinatorial optimization, such as a production scheduling problem, a vehicle routing problem, and a distributed network resource allocation optimization problem, etc.