A Multi-Objective Evolutionary Algorithm Model for Product Form Design Based on Improved SPEA2

: As a Kansei engineering design expert system, the product form design multi-objective evolutionary algorithm model (PFDMOEAM) contains various methods. Among them, the multi-objective evolutionary algorithm (MOEA) is the key to determine the performance of the model. Due to the deﬁciency of MOEA, the traditional PFDMOEAM has limited innovation and application value for designers. In this paper, we propose a novel PFDMOEAM with an improved strength Pareto evolutionary algorithm 2 (ISPEA2) as the core and combining the elliptic Fourier analysis (EFA) and the entropy weight and technique for order preference by similarity to ideal solution (entropy-TOPSIS) methods. Based on the improvement of the original operators in SPEA2 and the introduction of a new operator, ISPEA2 outperforms SPEA2 in convergence and diversity simultaneously. The proposed model takes full advantage of this superiority, and further combines the EFA method’s high accuracy and degree of multi-method integration, as well as the entropy-TOPSIS method’s good objectivity and operability, so it has excellent comprehensive performance and innovative application value. The feasibility and e ﬀ ectiveness of the model are veriﬁed by a case study of a car form design. The simulation system of the model is developed, and the simulation results demonstrate that the model can provide a universal and e ﬀ ective tool for designers to carry out multi-objective evolutionary design of product form.


Introduction
Rich and homogenous products mean manufacturers face increasingly serious competition. How to better meet consumers' increasingly diverse affective needs through product forms is the key to product design and the inevitable requirement for manufacturers to win this competition [1,2]. Kansei engineering or Kansei ergonomics has demonstrated an extremely strong relationship between customers' affective needs and product physical properties such as function and form [3]. As a Kansei engineering design expert system which can automatically generate product form according to consumers' diverse affective responses, the product form design multi-objective evolutionary algorithm model (PFDMOEAM) has increasingly significant value. A complete PFDMOEAM usually includes three modules: design analysis (DA), multi-objective optimization (MOO) and optimal solution selection (OSS). Each module contains a main method, i.e., product form and image analysis (PFIA) method, a multi-objective evolutionary algorithm (MOEA) and multi-attribute decision-making (MADM) method. Among all the methods, MOEA in the MOO module is the key to determine the comprehensive performance of the model, which mainly includes innovation, accuracy and operability. To improve the performance of PFDMOEAM, the traditional studies have tested different existing MOEAs, and can be divided into two kinds according to the MOEAs adopted. The first kind is to operable and can evaluate multiple objectives [16]. These characteristics are well in line with the objective of OSS module. Therefore, this paper employs entropy-TOPSIS as the MADM method in the OSS module.
The main contributions of this paper are summarized as follows: (1) To eliminate the major performance obstacle of PFDMOEAM, a novel ISPEA2 algorithm is first proposed based on the SPEA2 algorithm selected from renowned second-generation classical MOEAs. By improving the original operators in SPEA2 and introducing a new operator, ISPEA2 has better convergence and diversity than SPEA2, which lays a solid foundation for improving the performance of PFDMOEAM.
(2) Taking the ISPEA2 algorithm as the core and combining the EFA and entropy-TOPSIS method, we construct a novel PFDMOEAM. The model makes full use of ISPEA2 s algorithm advantages, and further combines the EFA method's high accuracy and degree of multi-method integration, as well as the entropy-TOPSIS method's good objectivity and operability, so it has excellent comprehensive performance in the multi-objective evolutionary design of product form.
(3) We applied the proposed PFDMOEAM to the car form design problem to demonstrate the innovation and application value of the model. The result proves that it can be used as an effective Kansei engineering design expert system and provides a useful tool for designers.
The rest of this paper is organized as follows. Section 2 presents a background overview of the related theory. Section 3 details the implementation procedures of the proposed model, and Section 4 demonstrates the discussion of the results. Finally, Section 5 provides some brief conclusions.

Research Framework
The overall framework of the proposed PFDMOEAM is presented in Figure 1, which consists of three modules, namely DA, MOO and OSS. The DA module is the first module of PFDMOEAM, which provides quantitative analysis data for the MOO module. The initial data contains two parts, product form data and product image data, which need to be obtained by using two sub-methods of the PFIA method, namely product form representation (PFR) and key image recognition (KIR). The PFR method primarily affects the innovation of the model by the accuracy and freedom of morphological representation, and the KIR method mainly affects the precision of the model through the accuracy and comprehensiveness of image recognition. Given the importance of data quality to the model, the PFIA method is often the basis and prerequisite for ensuring a good comprehensive performance of PFDMOEAM. Firstly, in the DA module of this paper, EFA is adopted as the PFIA method to get the principal component (PC) scores, and combined with image analysis, the image evaluation mean score (IEMS) of the target adjectives is obtained. Secondly, in the MOO module, the genetic algorithm and back propagation neural network (GABP) is first used to establish the non-linear mapping relationship between PC scores and IEMSs, namely multi-objective prediction model (MOPM), which is further used as the fitness function of ISPEA2. Then the Pareto solutions are derived through ISPEA2 improved by the improved crossover operator (ICO), the improved adaptive mutation operator (IAMO) and the modification operator (MO). In ISPEA2, ICO and IAMO are used to adjust the search space dynamically, and the newly introduced MO ensures the algorithm converges to global optimum efficiently. Through these operators, ISPEA2 successfully makes up for the defects of SPEA2 and retains its superiority of diversity. Finally, in the OSS module, the MADM method mainly completes accurate multi-objective evaluation of Pareto design solutions generated in the MOO module. The evaluation results can help designers to observe the actual effect of the model in time and determine the optimized direction of parameters in main methods. Due to its characteristics comply with the objectives of OSS module, entropy-TOPSIS is adopted as the MADM method to get the optimal solution.

Elliptic Fourier Analysis (EFA)
EFA [17] is a technique for converting point coordinate data into Fourier harmonic coefficients to quantify the profile of an object and has been widely used in many fields [18][19][20][21]. Taking into account all the information of the profile, EFA is better than other traditional morphometric approaches [22]. Its implementation flow is given below (Figure 2).

Calculation of Elliptic Fourier Descriptor (EFD)
After picture preprocessing, we get the closed product profile curve, which can be regarded as the motion trajectory of the moving point and expressed as two following functions of the arc length (l) from a point to the starting point on the profile curve [23]: where x(l) and y(l) are continuous piecewise linear functions; a 0 2 and c 0 2 are the abscissa and ordinate of the profile center point, respectively; n is the number of ellipse Fourier harmonics; N is the maximum number of ellipse Fourier harmonics; l is the arc length parameter; L is the profile perimeter; a n , b n , c n and d n represent the elliptic Fourier coefficients of the nth harmonic and can be defined as follows [17]: where K is the total number of sampling points of the profile; ∆x p and ∆y p are the displacements from the point p to the point p + 1 along the X-axis and Y-axis directions, respectively, then: Thus, product profile can be described by an elliptic Fourier coefficient matrix containing n harmonics, which is given as F is the EFD that approximates the profile curve by controlling the harmonic number n. When n takes different values, the fitting effect of profile curve will be different. In general, a few low-frequency harmonics maintain the global features of profile curve, and high-frequency harmonics affect its finer features.

Normalization of EFDs
Because the position, direction and scale of different product profiles are different after picture preprocessing, their EFDs are inconsistent and needs to be normalized to get the normalized elliptic Fourier descriptor (NEFD) F * [23]. Let the normalized coefficients of nth harmonic are a * n , b * n , c * n and d * n , then: (1) Since the DC component of EFD only reflects the position of the profile, it is not useful in describing the form. Let the DC component a * 0 = c * 0 = 0. The center of the first ellipse of the profile can be moved to the coordinate origin to achieve the normalization of the position.
(2) Transform the elliptic Fourier coefficients according to the following equations to achieve the normalization of direction: sin φ 1 cos φ 1 a n c n b n d n cos nθ 1 sin nθ 1 − sin nθ 1 cos nθ 1 (5) where θ 1 is the angle between the starting point of the profile and the main axis of the ellipse; ϕ 1 is the sequence of rotated ellipse; a 1 , b 1 , c 1 and d 1 are the four coefficients of the first ellipse.
(3) Calculate the size of the first ellipse E * 1 and divide each descriptor by E * 1 to achieve the normalization of scale:

Principal Component Analysis (PCA)
Principal component analysis (PCA) is an indirect ordination technique for obtaining a low dimensional representation of multivariate data [24]. We can get the PC scores of all product profiles by PCA. Using the PC scores as ordinary quantitative characters, the form can be quantitatively analyzed [17]. PC scores constitute crucial data for supporting product profile reconstruction, image analysis and multi-objective evolutionary calculation.

Image Analysis
As a quantitative analysis technique to study the affective needs of customers, the analysis of Kansei vocabularies is generally used to recognize the key image words influencing consumers' perception of product form [25,26]. This paper adopts this technique to recognize key image words. The steps are as follows: Step 1: Collect image words, that is, adjectives, and organize the subjects to conduct an informal survey to pick out appropriate adjectives.
Step 2: Invite experts to set up a focus team and use Delphi method to select the most representative adjectives.
Step 3: Combine the representative adjectives with the product samples to design a questionnaire. Conduct questionnaire survey to assess the experts' perception of the product image using the semantic differential scale and 7-point Likert scale which are two commonly used scales to quantify human perceptual interpretations [27]. All experts' evaluation scores for each product sample are averaged to achieve the IEMSs matrix.
Step 4: After concatenating the IEMSs matrix and the PC scores matrix to create a correlation analysis matrix, correlation analysis is performed on the matrix to obtain correlation coefficients. The correlation coefficient is then brought into the image objective comprehensive scoring model (IOCSM) to calculate the final image objective comprehensive score (FIOCS). For more detailed information about the IOCSM, please see Wang et al. [14].
Step 5: Determine the target adjectives and their rankings according to the FIOCS obtained in Step 4.

Genetic Algorithm and Back Propagation (GABP)
In the MOPM of MOO module, the PC scores and IEMSs, which are input and output variables, respectively, are continuous data. Wang et al. [28] demonstrated that the prediction and fitting ability of neural network model for continuous data is more stable and stronger than the multiple linear regression model used in current authors' previous research [14], and the combination of GABP can further better predict the nonlinear mapping relationship between input and output variables [6,29,30].
Consequently, this paper uses GABP as the prediction technique in MOPM to achieve the high-precision prediction requirements of MOEA. For detailed implementation procedures of the technique, please see, for example, Huang et al. [30].

Improved Strength Pareto Evolutionary Algorithm 2 (ISPEA2)
The classical SPEA2 algorithm uses the nearest-neighbor density estimation technique that preserve the diversity of non-inferior solutions, the hybrid selection technique that maintains the evolutionary direction, and the environmental selection technique that retains the Pareto solutions [31]. The proposed ISPEA2 keeps these characteristics of SPEA2, and obtains better convergence and diversity by the ICO, IAMO and MO. Its implementation flow is demonstrated in Figure 3, which consists of seven main steps: Step 1: Initialize population P 0 and empty archive (external set) A 0 = ϕ. Set the generation number t = 0.
Step 2: The neural networks trained by GABP are used to calculate the fitness values of the multi-objective images corresponding to them. The fitness values are then assigned to each individual in P t and A t . Specifically, each individual in P t and A t is first given an intensity value St(i) according to dominance relation, where St(i) can be defined as: where the symbol represents the Pareto dominance relation. Then, calculate the raw fitness based on each individual's intensity value, which is given below: Then, calculate the density of the individual to distinguish individuals with the same raw fitness: where σ γ i denotes the distance from individual i to the γth adjacent individual after an ascending order of distance between individual i and all other individuals. Finally, the fitness of individual i can be expressed as follows: Step 3: Identify all non-dominant individuals and copy them from P t and A t to A t+1 using an environmental selection technique. If the size of A t+1 exceeds its maximum permitted size, fill A t+1 with the best non-inferior solutions by the truncation operation. Otherwise, fill A t+1 with the dominant individuals in P t and A t .
Step 4: If the stopping criterion is satisfied then use the non-dominated individuals in A t+1 as the final Pareto solution set and terminate.
Step 5: Fill the mating pool with a binary tournament selection to replace A t+1 .
Step 6: Apply the ICO, IAMO and MO to obtain the resulting population P t+1 . The calculation process of each operator is as follows: (1) Let PC scores vectors g t i and k t i be two parent individuals of the ith generation to be genetically manipulated after the operation of the above steps, and two descendant individuals after the crossover operation through ICO are expressed as where α is the random number with uniform distribution on [−0.5, 1.5], and α + β = 1; ξ is the number of PC.
(2) Let PC scores vector h t i be the parent individual of the tth generation, the descendant individual after the mutation operation through IAMO is expressed as where ms is the coefficient of mutation scale; md is the decay rate of mutation; cn is the current number of iteration; mn is the maximum number of iterations; rn 1 is a random number of standard normal distribution; ub and lb are the upper and lower limits of the boundary of a gene, respectively.
(3) The MO modifies the results of the crossover and mutation operations with a certain modification probability according to the following equation: where f t+1 ij represents the descendant individuals obtained after the crossover and mutation operations, i is a set of all PC numbers; i * corresponds to the PC numbers with consistent correlation coefficient between PC scores and IEMSs of target adjectives; mc i * is the modification coefficient of the i * th PC number; rn 2 is a random number in the interval [4,5]; λ i * z is the correlation coefficient of the zth target adjective of the i * th PC number, and it is considered that there is a correlation between two variables when [32]; Z is the number of target adjectives. By replacing the gene corresponding to the position of i * of f t+1 ij with f t+1 i * j , the modification operation can be completed.
Each operator acts on the algorithm as follows. Firstly, the value range of α and β in ISPEA2 is not limited to [0,1] in SPEA2, this ensures that the search space of the ICO in ISPEA2 covers all the neighborhoods of g t i and k t i . Therefore, ICO can make the probability of searching for more different offspring higher and the diversity of ISPEA2 better. Secondly, the IAMO ensures that the mutation rate changes dynamically with the increase of the number of iterations. In the initial stage of the iteration, the global optimal solutions are searched by a large mutation rate, and at the end of the iteration, local optimal solutions are fully searched by a small mutation rate. This improves the convergence of ISPEA2. Finally, the MO helps to achieve the goal of continuously improving the scores of target adjectives in iterations. Since the value range of λ i * z and rn 2 determines the reasonable range of mc i * as the modification coefficient, which ensures that no matter positive and negative correlation, individual genes consisting of PC scores of the position of the i * th PC number can be evolved towards the direction of increasing the scores of target adjectives with the pre-set modification probability in each iteration. Therefore, the introduction of MO can significantly improve the convergence of the algorithm. By reasonably presetting the modification probability, it can effectively maintain a good balance between the convergence and diversity of SPEA2.
Step 7: Increase the generation number t by 1 and go to Step 2.

Entropy Weight and Technique for Order Preference by Similarity to Ideal Solution (Entropy-TOPSIS)
Entropy-TOPSIS determines the weight of each target adjective based on the entropy weight method, and evaluates Pareto solutions and gets the optimal solution based on the TOPSIS technique. By calculating the comprehensive evaluation index Ω i , all Pareto solutions can be ranked in order. For further details and analysis, please see, for example, Chauhan et al. [33].

Empirical Study
The feasibility and effectiveness of the proposed model is demonstrated by a case study of car form design. The detailed implementation procedures are described in the following sections.

Design Analysis
Cars are mature and diverse products that are familiar to consumers. They are ideal research objects in the field of product form design [34,35]. As the most crucial form feature of a car, a car profile is mainly represented by a closed curve including the engine line, the windshield line, the roof line, the wheelbase and the wheel arch line [36]. In view of the fact that most researchers have used a 2D profile to define the form feature of a product [37], this study verifies the proposed model by taking a 2D car profile as an example. In consideration of the familiarity and satisfaction of customers, we first collected pictures of 125 mid-size and three-box cars from 32 different brands in the Chinese market. Subsequently, a focus group of nine experts (five men and four women) with more than two years of car design experience were formed to classify and select representative cars based on their similarity in visual features. By using the Kawakida Jirou (KJ) method, multidimensional scaling analysis (MDS) and cluster analysis [38,39], 45 representative samples were selected to form the final benchmark car sample set. Figure 4 demonstrates the sample pictures, and Table 1 manifests their brands and models.  After the image preprocessing process shown in Figure 5, we got the data sequence of point coordinates of representative samples' profiles.
As can be seen from Figure 6, the tire sizes of different representative samples have not changed significantly, so their influence on the product image is negligible. Consequently, the extracted car profiles do not contain tire information (see Figure 7).
We performed EFA on data sequence of point coordinates of the profiles according to the Equations (1)- (9) and obtained the F * of each sample. Table 2 demonstrates the F * for sample No.1.
To observe the effect of profile reconstruction, sample profiles were reconstructed by inverse Fourier transform of F * (see Figure 8).
By comparison with the original profile, n was determined to be 64, because at this harmonic number, the reconstructed profile had most of the features of the original profile and the noise was very small. Then, the F * of each sample was reshaped from a matrix of dimensions 64 × 4 into a row vector, and the matrices of all samples were assembled into a matrix of dimensions 45 × 256. After that, we performed PCA on the matrix to obtain a low-dimensional PC scores matrix.
To ensure high-accuracy reconstruction of the sample profile and reduce the computational complexity of the proposed model, the number of PC was set to 7 based on the cumulative variance contribution rate. At this point, the cumulative variance contribution rate exceeded 99.95% (see Table 3) and the PC scores matrix was a matrix of dimensions 45 × 7 ( Table 4, whole set of data see Table S2).

Image Analysis
Firstly, 128 affective adjectives describing car profiles were collected from literature, car catalogs and professionals. Next, an informal survey was conducted on 16 industrial design students (8 males and 8 females) to select 36 relevant adjectives preliminarily (see Table 5), which were suitable for these subjects to express their perceptions on the representative sample profiles. Then, the focus group finalized 11 most representative adjectives, with application of the Delphi method. After that, questionnaires were designed by combining 45 representative samples with the 11 representative adjectives, using a 7-point semantic differential method (-3-3), in which −3 and 3 are the lowest and the highest score of adjectives, respectively. Subsequently, the IEMSs was calculated after the focus group completed the questionnaire ( Table 6, whole set of data see Table S3), and the correlation analysis matrix for all samples was obtained by concatenating the IEMSs and PC scores matrixes. Finally, the FIOCS and rankings of 11 representative adjectives were obtained via using the IOCSM ( Table 7). The top three adjectives in the rankings were set as the target adjectives in the case study of this paper, namely "luxurious", "steady", and "modern".

Construction of Multi-Objective Prediction Model (MOPM)
In this section, we constructed three BPNNs optimized by GA to predict calculated scores of three target adjectives. Each BPNN selected the classic three-layer structure, and the number of neurons in the input layer, hidden layer and output layer were 7, 10 and 1, respectively. GA was adopted to generate the optimal initial weights and thresholds for three BPNNs. After being trained, these BPNNs were further used as the fitness function to predict the calculated scores of three target adjectives in ISPEA2. Table 8 reflects the specific parameter values of GABP, some of which are derived from Zhang et al. [40]. The details of the algorithm parameters are manifested in Table 8, The values of some parameters in Table 8 and Table 11 are referenced to the studies. Figure 9 demonstrates the prediction errors of three BPNNs when the number of iterations is 20,000.
As can be seen from Figure 9, there are small errors between the predicted and true values of all samples for three target adjectives. In addition, we further used the root mean square error (RMSE) to evaluate the performances of these BPNNs, which is defined as where ns is the number of samples; ψ * i and ψ i are the BPNN predicted and true values of the ith sample, respectively. Table 9 reflects that RMSE of the three BPNNs decrease rapidly as the number of iterations increases. After 20,000 epochs, the final RMSE values are 0.0455, 0.0534 and 0.0526, respectively, all within the acceptable tolerance of less than 0.1, which demonstrates that the convergence performance of each network is satisfactory. It can be seen that GABP has built a MOPM with high prediction accuracy, which provides a strong guarantee for the effective implementation of ISPEA2.

Pareto Solutions Generation
In this section, we first determined the key parameters of ISPEA2. Table 10 illustrates the correlation coefficients between PC scores and IEMSs of the three target adjectives, which indicate that PC1 and PC2 are consistent positive and negative correlation with the three target adjectives, respectively. According to the Equation (19), m 1 and m 2 were calculated as 0.718 and −0.54 by using the correlation coefficients of PC1 and PC2 in Table 10. Then the modification coefficients mc 1 and mc 2 were determined by combining Equation (18) and random number rn 2 . The purpose of setting the modification probability is to accelerate the convergence of ISPEA2 without affecting the diversity of the algorithm too much, so as to ensure a good balance between the convergence and diversity. As such, it is necessary to compare the scores of multi-object image of Pareto solutions generated after an independent run of the algorithm and the visual difference of the profile shape of the generated solutions. The higher the scores and the larger the shape difference, the better the parameter value will be. Through multiple simulation experiments on different modification probability values, we observed the scores and visual differences of the Pareto solutions generated in different independent runs, and finally set the modification probability to 0.3. The parameter values of ISPEA2 are given in Table 11, some of which are derived from Wen et al. [41].  After running ISPEA2 according to the above parameters, we obtained a total of 30 Pareto solutions (design alternatives), whose PC scores and target adjective calculation scores are shown in Table 12 (whole set of data see Table S4).

Optimal Solution Selection and Verification
In order to obtain the optimal solution, we calculated the weights of three target adjectives and the comprehensive evaluation index Ω i of all Pareto solutions by using entropy-TOPSIS, Table 13 manifests the weights of three target adjectives. Table 14 demonstrates the detailed Ω values and rankings of 30 Pareto solutions. Since Ω 7 ranks first, it was determined to be the optimal solution. To verify the result, four representative Pareto solutions in the Pareto solution set were selected at the starting point, the middle adjacent two points, and the ending point, respectively, which are marked out with red solid dots in Figure 10. The PC scores, target adjective calculation scores, Ω i values and profiles of the four selected Pareto solutions are demonstrated in Table 15. After combining the profiles of these solutions with three target adjectives to create a questionnaire, the focus group was invited again to evaluate these solutions. The result is manifested in Table 16, from which we can see that the expert evaluation result is consistent with the automatic calculation result of the model in this paper.

Discussion
In this paper, we have proposed a PFDMOEAM which integrates three main methods: ISPEA2, EFA and entropy-TOPSIS. The comprehensive performance superiorities of the model derive from the performance benefits of these methods.

Performance of ISPEA2
The performance of MOEA algorithm is the most critical factor in determining the comprehensive performance of PFDMOEAM. After the comparison of second generation classical MOEAs, we chose to improve SPEA2 into the MOEA required by the PFDMOEAM built in this paper. To verify the performance of ISPEA2, ISPEA2 and SPEA2 were independently run 30 times under the same parameter settings (see Table 11) to compare the performance of two algorithms by using convergence metric, spacing metric and coverage metric [42][43][44]. Since the true Pareto front in this study is unknown, we adopted a common comparative method suitable for this situation, which is to construct a reference Pareto front to represent the true Pareto front by a pair-wise comparison after collecting all Pareto solutions from all independent runs of the two algorithms [7,45]. Let A = (A 1 , A 2 , . . . , A |A| ) be a solution set obtained by ISPEA2 or SPEA2, and P * = (P * 1 , P * 2 , . . . , P * |P * | ) be a solution set on the true Pareto front, the three metrics can be calculated as follows.
Convergence metric C represents the average distance from the solution set obtained to the true Pareto front, which can be calculated as [42,46,47]: where d i is the minimum normalized Euclidean distance from each solution in A to P * ; f max z and f min z are the maximum and minimum function values of the zth objective function in the true Pareto solution set P * , respectively. The smaller the C value is, the better the convergence performance will be. When C = 0, it means that all solutions in A are on the true Pareto front. Spacing metric S denotes the distance between two consecutive solutions and measures the diversity performance of algorithm, which is computed as [43,47]: where d is the mean of d i . A smaller S value indicates that the obtained solutions are distributed more evenly throughout the objective space.
Coverage metric C(A 1 , A 2 ) denotes that the obtained solutions of one algorithm dominate the obtained solutions of the other algorithm, which is defined as [44,48]: where C(A 1 , A 2 ) means the ratio of the number of solutions dominated by A 1 in A 2 to all solutions in A 2 , and C(A 2 , A 1 ) indicates the opposite case. If C(A 1 , A 2 ) ≥ C(A 2 , A 1 ), then A 1 are better than A 2 . To compare C(A 1 , A 2 ) values of the two algorithms from a statistical perspective, A 1 and A 2 solution sets are obtained by further calculation of Pareto solution generation after collecting all solutions of all runs of the two algorithms, respectively. After calculating the metrics mentioned above, we computed the means and standard deviations of C and S metrics and the C(A 1 , A 2 ) metric values of two algorithms (Table 17) and plotted the trend graph of C and S metric values of 30 independent runs ( Figure 11).  From Table 17 and Figure 11, we can see that the C metric value of ISPEA2 is obviously smaller than that of SPEA2 and the S metric value of ISPEA2 is also lower than that of SPEA2, which manifests that the convergence of ISPEA2 has been obviously improved, and the diversity of ISPEA2 has also been enhanced. In addition, C(A 2 , A 1 ) = 0 means that there is no solution in SPEA2 that can dominate the solution in ISPEA2. The high value of the C(A 1 , A 2 ) metric demonstrates that the outcomes of ISPEA2 remarkably dominate the outcomes of SPEA2, and the solutions of ISPEA2 are closer to the true Pareto front and more widely distributed. Consequently, the comparative results prove that ISPEA2 has better convergence and diversity. It should be pointed out that the standard deviations of two metrics are very close, which demonstrates that the stability of two algorithms is relatively consistent.
The ISPEA2 algorithm is the main source of the proposed PFDMOEAM's comprehensive performance superiority and value, helping the model to generate high quality Pareto design solutions. The overall forms of these solutions are rich and free, and the form details are fine and smooth, even the morphological information of local structures such as car lights can be observed. Rich and fine forms match the multi-image objective better, which significantly improves the innovation and application value of the proposed PFDMOEAM.

Superiorities and Values of EFA for Product Form Design Multi-Objective Evolutionary Algorithm Model (PFDMOEAM)
In the authors' previous work, it has been proved that EFA as an effective PFIA method has shown significant superiorities in the DA module. Compared with other traditional PFIA methods, it improves the precision and degree of freedom of PFR, and enhances the accuracy and comprehensiveness of KIR. This paper further proves that EFA can provide a new idea and way for optimizing the performance of MOEA in the MOO module. This is reflected in the fact that the MO operator is proposed by using the consistency correlation between PC scores and IEMSs both obtained by EFA. Because PC scores can reflect overall form features, and the perception mode of product image also adheres to the principle of overall form priority of product image recognition [49], the relationships of PC scores and IEMSs relative to product form are both more based on the whole rather than the local. By using these relations, the MO operator adjusts the evolutionary direction, which reduces the amount of computation and improves the efficiency and performance of ISPEA2 algorithm in approximating the true Pareto solution set. It can be seen that EFA not only promotes the achievement of the goals of DA and MOO modules better, but also creates an organic connection between different methods in these two modules, which helps to integrate the model as a whole more effectively. Therefore, EFA is a significant basis for gaining the superiorities of the proposed PFDMOEAM.

Complete PFDMOEAM and Entropy-TOPSIS
As the last module of the complete PFDMOEAM, the main task of the OSS module is to use the MADM method to evaluate and rank Pareto design solutions according to multi-objective calculated scores. Based on the superiorities of innovation and accuracy gained in the DA and MOO modules, the goal of OSS module is to consolidate these superiorities and get accurate MADM results. Therefore, an objective, accurate and operable MADM method is needed in this module. Compared with other subjective MADM methods, the entropy-TOPSIS method has better objectivity, accuracy and operability, so it can achieve the goal of the OSS module more effectively. The consistency of the results of Tables 15 and 16 demonstrates that the entropy-TOPSIS method has good applicability to PFDMOEAM and has been integrated with the other main methods to provide strong support for the comprehensive performance benefits of the model. Based on the above methods, the automatic simulation system of the model was developed by using MATLAB software. The system provides a visual interface that allows designers to carry on designing and observing the results intuitively and efficiently ( Figure 12). The implementation and operation of simulation system proves that the proposed model can be used as a universal Kansei engineering design expert system with great innovation and application value.
Note that 45 representative samples were used in the case study for test, but the proposed PFDMOEAM can increase or decrease the sample size according to actual needs. Since the PCA in the EFA method requires a certain number of samples, the sample size should not be too small. In theory, the larger the sample size, the more reliable the result will be.

Conclusions
To overcome the major obstacle that restricts the performance of traditional PFDMOEAMs, we have proposed a novel ISPEA2 algorithm. Taking ISPEA2 as the core and combining the EFA and entropy-TOPSIS methods, we have further constructed a novel PFDMOEAM and developed its simulation system in this paper. The simulation results verify the feasibility and effectiveness of the model.
(1) Since MOEA is the key to determining the performance of PFDMOEAM, we have proposed a novel ISPEA2 algorithm based on the SPEA2 algorithm selected from the renowned second generation classical MOEAs. In addition to improving the original operators in SPEA2, the MO operator was introduced for the first time to optimize the performance of ISPEA2. Through these operators, ISPEA2 is superior to SPEA2 in convergence and diversity, which lays a solid foundation for improving the comprehensive performance of the model.
(2) For a better complete model, we have further constructed a novel PFDMOEAM that combines ISPEA2, EFA and entropy-TOPSIS. It makes full use of the superiorities of these methods to maximize its comprehensive performance, which are the excellent convergence and diversity of ISPEA2 as MOEA algorithm, the high accuracy and degree of multi-method integration of EFA as PFIA method, and the good objectivity and operability of entropy-TOPSIS as the MADM method. Therefore, the model exhibits excellent comprehensive performance in multi-objective evolutionary design of product form.
(3) A case study of car form design has been used to elaborate the process of the proposed model, and its simulation system has been developed. The simulation results demonstrate that the model can generate a variety of new car profiles that meet the multi-image objective, offer powerful support for the scientific MADM of car form optimization, and provide a universal and effective design tool for designers. Although this study takes car form as an example, the model is also applicable to the form design of other products.
It is acknowledged that the proposed model needs to integrate many different methods, which makes the debugging process a bit complicated. However, this problem has been properly implemented in practice. In addition, this model is only used in the design of a 2D product profile. Further exploring the application of the model in 3D product form design and combining the model with product color design will be the focus of future research.