Shape Design Optimization of a Robot Arm Using a Surrogate-Based Evolutionary Approach

: In the design optimization of robot arms, the use of simulation technologies for modeling and optimizing the objective functions is still challenging. The di ﬃ culty is not only associated with the large computational cost of high-ﬁdelity structural simulations but also linked to the reasonable compromise between the multiple conﬂicting objectives of robot arms. In this paper we propose a surrogate-based evolutionary optimization (SBEO) method via a global optimization approach, which incorporates the response surface method (RSM) and multi-objective evolutionary algorithm by decomposition (the di ﬀ erential evolution (DE ) variant) (MOEA / D-DE) to tackle the shape design optimization problem of robot arms for achieving high speed performance. The computer-aided engineering (CAE) tools such as CAE solvers, computer-aided design (CAD) Inventor, and ﬁnite element method (FEM) ANSYS are ﬁrst used to produce the design and assess the performance of the robot arm. The surrogate model constructed on the basis of Box–Behnken design is then used in the MOEA / D-DE, which includes the process of selection, recombination, and mutation, to optimize the robot arm. The performance of the optimized robot arm is compared with the baseline one to validate the correctness and e ﬀ ectiveness of the proposed method. The results obtained for the adopted example show that the proposed method can not only signiﬁcantly improve the robot arm performance and save computational cost but may also be deployed to solve other complex design optimization problems.


Introduction
Nowadays, facing the megatrend for smart factories of the future, industrial robots play a greatly important role not only in traditional pick-and-place tasks but also in many of the precision applications such as assembly, welding, and machining. These robots, in order to satisfy the requirements of new applications, must be designed to meet high end/performance requirements, i.e., high accuracy and high efficiency. The design process of an industrial robot arm with high speed performance is still challenging due to the interconnecting roles of adjacent joints, which may significantly impact the efficiency and stability of the robot arm. For example, the deformation and residual vibration of the robot arm will cause its positioning accuracy to decrease, and the robot arm mass and moment of inertia will affect its dynamic behavior, such as controllability and efficiency; moreover, these factors generally are in conflict with each other. Therefore, design optimization of robot arms is becoming increasingly important.
With the advance of digital twins in smart factories, the use of simulation technologies to address computation-expensive and multi-objective optimization issues has become an important research topic. Design optimization is the process of establishing the right combination of product parameters that satisfy project requirements [1]. Structural optimization can be separated into three main categories: size, shape, and topology optimization [2]. In a structural design process, a topology optimization is first performed to obtain an initial sketch of the optimal configuration of the structure; the suggested configuration is interpreted to draw an engineering design, and this design is then optimized using detailed sizing and shape optimization methods together with realistic design constraints. These methods are not only able to optimize the robot arm but can also be extended to the design of the whole robot system. For the purpose of improving the performance of industrial robots, a number of research works have been carried out with respect to the related topics in recent years. Some works focused on the selection criteria of the servo motor and gearbox for optimizing the servo drive system [3,4]. A few reports studied new control algorithms to obtain high speed performance, such as residual vibration suppression [5]. Certain researchers used the finite element method (FEM) to optimize robot structure. Roy et al. used the finite element method to design a high-performance robot arm and reported experimental structural data to support the design of the robot arm [6]. Park and Asada studied the design optimization of the mechanical structure and control of a two-link robot for which the weight and inertia of the robot arm were minimized to attain high-speed positioning [7]. Sahu and Choudhury used the FEM of the ANSYS workbench to optimize the performance of a six-axis industrial robot via stress, modal, and dynamic behavior analyses [8][9][10]. Pupăză et al. also used an FEM-based optimization technique to study the influence of the static and dynamic behaviors on the positioning accuracy of robots [11]. Chen et al. analyzed the static strength and stiffness of a two-link robot and then used a topology optimization method to design the structure for quality improvement [12]. In the above studies, it was revealed that CAE together with FEM-based design optimization can be an important tool for attaining structural performance improvement in real-world engineering applications. However, for cases with large design space or the requirement of attaining numerous design samples, the high computational effort of analysis may hinder the implementation of the conventional FEM-based design optimization methods in the design optimization process. On the other hand, only a few researchers have studied the multi-objective and computational inefficiency problems related to the design optimization of industrial robots.
The response surface method (RSM), also known as the surrogate model or approximation model, can be used to improve the design efficiency. In RSM, a set of sample data (obtained by running the computer code/simulation) are used to build the surrogate model, which are sufficient to predict the output at untried points within the design space. Because the predictions obtained from a surrogate model are generally much more efficient than those from numerical analysis code, the computational cost associated with each search in RSM is therefore negligible. Surrogate-model-based optimizations are of particular interest for engineering design when coupled with high-fidelity and expensive analysis codes, such as computation fluid dynamics (CFD) and computational structural dynamics (CSD) [13]. In many studies, RSM together with design of experiment (DOE) was used to construct the surrogate models and find the best-fitting solutions [14][15][16][17]. Some researchers applied the adaptive RSM to search for the global optima of complex mechanical systems [18][19][20]. However, in these cases, they did not consider the situation of multi-objective functions with conflicting properties [21] which may lead to being trapped into a local optimum solution during the design optimization process.
In this paper, the shape design optimization of robot arms for high speed performance via the surrogate-based evolutionary approach is investigated. First, based on the CAD model of a robot arm, the design variables, constraints, and objective functions of the optimization problem are defined. Next, a quadratic RSM together with Box-Behnken design (BBD) and FEM is used to construct the surrogate model. In solving the optimization problem, the multi-objective evolutionary algorithm by decomposition (DE variant) (MOEA/D-DE) technique which includes evolutionary operators of selection, recombination, and mutation is employed to address the multiple objective functions with conflicting features and search for the Pareto-optimal sets of global optimization. The weighted sum method is finally used to determine the best candidate among the generated population samples through fitness checking. The proposed paper is able to make two significant contributions: one is to provide a surrogate-based evolutionary optimization approach to deploy in the robot arm design optimization process to obtain higher performance and efficiency, as well as to reduce computational cost; the other one is to solve the shape design optimization problem of robot arms considering multiple, conflicting objectives for higher accuracy and controllability simultaneously.
The rest of the paper is organized as follows. The techniques for robot arm analysis and design, which involve the target performance to be achieved and the framework of the proposed surrogate-based evolutionary optimization (SBEO) method, are described in Section 2. The procedure for the optimal shape design of a specific robot arm is given in Section 3. Some optimal-design-related simulation results, such as Pareto front history, improvement of performance, computational efficiency, and sensitivity analysis, are given in Section 4. Section 5 presents a discussion on the advantages, limitations, and future works of the proposed method.

Architecture of the Industrial Robot
The PMC6VA030 robot, the six-axis vertical articulated industrial manipulator shown in Figure 1, is fabricated by PMC, Taiwan. The maximum payload and reach of the robot are 30 kg and 1800 mm, respectively. The architecture of the robot is composed of six revolute joints and six robot arms. Each robot arm is driven by a joint module which consists of a servo motor, reduction gear, and encoder.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 17 conflicting features and search for the Pareto-optimal sets of global optimization. The weighted sum method is finally used to determine the best candidate among the generated population samples through fitness checking. The proposed paper is able to make two significant contributions: one is to provide a surrogate-based evolutionary optimization approach to deploy in the robot arm design optimization process to obtain higher performance and efficiency, as well as to reduce computational cost; the other one is to solve the shape design optimization problem of robot arms considering multiple, conflicting objectives for higher accuracy and controllability simultaneously. The rest of the paper is organized as follows. The techniques for robot arm analysis and design, which involve the target performance to be achieved and the framework of the proposed surrogatebased evolutionary optimization (SBEO) method, are described in Section 2. The procedure for the optimal shape design of a specific robot arm is given in Section 3. Some optimal-design-related simulation results, such as Pareto front history, improvement of performance, computational efficiency, and sensitivity analysis, are given in Section 4. Section 5 presents a discussion on the advantages, limitations, and future works of the proposed method.

Architecture of the Industrial Robot
The PMC6VA030 robot, the six-axis vertical articulated industrial manipulator shown in Figure  1, is fabricated by PMC, Taiwan. The maximum payload and reach of the robot are 30 kg and 1800 mm, respectively. The architecture of the robot is composed of six revolute joints and six robot arms. Each robot arm is driven by a joint module which consists of a servo motor, reduction gear, and encoder.

Problem Description of the Robot Arm
The accuracy, controllability, and efficiency are important performance aspects of a high-speed industrial robot. Due to the robot arms playing roles of interconnection between the adjacent joints, their characteristics are able to influence the system performance very well. When the power capacities of the joint modules are fixed, the shape dimensions of the robot arms are able to significantly affect the robot performance. Therefore, shape design optimization of the robot arms to obtain the minimum weight, maximum stiffness, and minimum deformation must be accomplished [22]. Figure 2a is a detailed engineering drawing of the baseline 2nd robot arm which includes all of the shape dimensions in it. A schematic view of the proposed 2nd robot arm with the chosen design variables is shown in Figure 2b. The unspecified dimensions on this sketch are all the same as those in Figure 2a and will not be optimized in this article under real engineering consideration.

Problem Description of the Robot Arm
The accuracy, controllability, and efficiency are important performance aspects of a high-speed industrial robot. Due to the robot arms playing roles of interconnection between the adjacent joints, their characteristics are able to influence the system performance very well. When the power capacities of the joint modules are fixed, the shape dimensions of the robot arms are able to significantly affect the robot performance. Therefore, shape design optimization of the robot arms to obtain the minimum weight, maximum stiffness, and minimum deformation must be accomplished [22]. Figure 2a is a detailed engineering drawing of the baseline 2nd robot arm which includes all of the shape dimensions in it. A schematic view of the proposed 2nd robot arm with the chosen design variables is shown in Figure 2b. The unspecified dimensions on this sketch are all the same as those in Figure 2a and will not be optimized in this article under real engineering consideration.

Objective Functions
In this paper, the weight (W), moment of inertia (J), and deformation (D) of the 2nd robot arm were chosen as tri-objective functions. The mathematical model of the multi-objective optimization problem is as follows [23]: where F(X), h(X), and g(X) are, respectively, the objective functions and the equality and inequality constraint functions; X ∈ R n is a vector of n design variables; p and q are the numbers of constraint functions. The multi-objective optimization problem can be expressed as a single equation via weighted sum and normalization: where j = 1, 2, 3, …., ps, and ps is the population size. w1, w2, and w3 are weighted factors for preference criteria.

Design Variables
Several independent geometric parameters of the 2nd robot arm, including shape dimension initial point Hi, position of inflection point Hm and Vm, radius of opening Rh, and arm thickness Ta, were chosen as the design variables. The ranges of the design variables in the evolutionary models are shown in Table 1.

Objective Functions
In this paper, the weight (W), moment of inertia (J), and deformation (D) of the 2nd robot arm were chosen as tri-objective functions. The mathematical model of the multi-objective optimization problem is as follows [23]: where F(X), h(X), and g(X) are, respectively, the objective functions and the equality and inequality constraint functions; X ∈ R n is a vector of n design variables; p and q are the numbers of constraint functions. The multi-objective optimization problem can be expressed as a single equation via weighted sum and normalization: where j = 1, 2, 3, . . . ., ps, and ps is the population size. w 1 , w 2 , and w 3 are weighted factors for preference criteria.

Design Variables
Several independent geometric parameters of the 2nd robot arm, including shape dimension initial point H i , position of inflection point H m and V m , radius of opening R h , and arm thickness T a , were chosen as the design variables. The ranges of the design variables in the evolutionary models are shown in Table 1. For statistical calculations, X i is related to x i via the following equation: where X i is the coded value of a dimensionless design variable which can be easily calculated; x i is the actual value of a design variable; x 0 is the actual value of a design variable at the center point; and x iU and x iL are the upper and lower bounds of a design variable.

Design Constraints
The above objective functions are subject to design constraints which are related to the geometry and mechanical behavior of the robot arm.

•
The first set of constraints is the limitation of the feasible design region for design variables based on Table 1.
• The second set of constraints considers the relationship between the radius R c of the circumcircle and design variables H i , H m , and V m . According to the formulation of the circumcircle, the coordinates (x, y) of the circumcenter of a triangle and its radius R c can be obtained as follows.
x = Hence, the constraint involving y should be • One constraint is related to the weight-to-stiffness ratio (W2SR), which is a measure of material mass efficiency used to provide unit resistance to elastic deformation. An improved W2SR enables the realization of longer spans without incurring a significant cost penalty, which is expressed in terms of material usage and deflections due to weight [24]. The W2SR computed from the weight and stiffness of a robot arm can be expressed as where F y is the applied force at the distal-end joint.

•
Another constraint is related to the weight-to-moment-of-inertia ratio (W2JR), a measure of a structural property. A higher W2JR will push resonance and anti-resonance peaks to lower frequencies, shrinking the operating bandwidth of the machine. The measure of W2JR in this paper is obtained as The above two measures of W2SR and W2JR are both "the smaller the better" when the design optimization of robot arms is targeting higher accuracy and controllability. Two additional constraints imposed on the objective functions are as follows:

Performance Analysis Using a CAE Solver
The performance of the robot arm was analyzed using a CAE solver in this paper. The model of an industrial robot built using Autodesk Inventor 3D CAD software was imported into the ANSYS MechanicalANSYS Parametric Design Language (APDL )FEM solver in the "*.STP" format. The material of the robot arm is aluminum 6061 with modulus of elasticity E = 69 GPa, density ρ = 2700 kg/m 3 , and Poisson's ratio ν = 0.3, applied force at the distal-end joint F Y = 1300 N, moment M X = 277 Nm, and M Z = 512 Nm.
The 2nd robot arm, together with the subsequent arms and distal-end joint, is shown in Figure 3a. All the subsequent arms are treated as point masses (also known as remote forces) which are used to transfer the force at the distal-end joint (right side) to the 2nd robot arm. Consequently, a resultant force F Y and two resultant moments M X and M Z are applied to the 2nd robot arm as shown in Figure 3b.

Performance Analysis Using a CAE Solver
The performance of the robot arm was analyzed using a CAE solver in this paper. The model of an industrial robot built using Autodesk Inventor 3D CAD software was imported into the ANSYS MechanicalANSYS Parametric Design Language (APDL )FEM solver in the "*.STP" format. The material of the robot arm is aluminum 6061 with modulus of elasticity E = 69 GPa, density ρ = 2700 kg/m 3 , and Poisson's ratio ν = 0.3, applied force at the distal-end joint FY = 1300 N, moment MX = 277 Nm, and MZ = 512 Nm.
The 2nd robot arm, together with the subsequent arms and distal-end joint, is shown in Figure  3a. All the subsequent arms are treated as point masses (also known as remote forces) which are used to transfer the force at the distal-end joint (right side) to the 2nd robot arm. Consequently, a resultant force FY and two resultant moments MX and MZ are applied to the 2nd robot arm as shown in Figure  3b. The FEM mesh of the 2nd robot arm is shown in Figure 4a, and the supporting boundary is placed at the proximal end (left side) of the robot arm with the harmonic reducer treated as a fixed connection as shown in Figure 4b. The FEA results are shown in Figure 5.   The FEM mesh of the 2nd robot arm is shown in Figure 4a, and the supporting boundary is placed at the proximal end (left side) of the robot arm with the harmonic reducer treated as a fixed connection as shown in Figure 4b. The FEA results are shown in Figure 5.

Performance Analysis Using a CAE Solver
The performance of the robot arm was analyzed using a CAE solver in this paper. The model of an industrial robot built using Autodesk Inventor 3D CAD software was imported into the ANSYS MechanicalANSYS Parametric Design Language (APDL )FEM solver in the "*.STP" format. The material of the robot arm is aluminum 6061 with modulus of elasticity E = 69 GPa, density ρ = 2700 kg/m 3 , and Poisson's ratio ν = 0.3, applied force at the distal-end joint FY = 1300 N, moment MX = 277 Nm, and MZ = 512 Nm.
The 2nd robot arm, together with the subsequent arms and distal-end joint, is shown in Figure  3a. All the subsequent arms are treated as point masses (also known as remote forces) which are used to transfer the force at the distal-end joint (right side) to the 2nd robot arm. Consequently, a resultant force FY and two resultant moments MX and MZ are applied to the 2nd robot arm as shown in Figure  3b.

Performance Analysis Using a CAE Solver
The performance of the robot arm was analyzed using a CAE solver in this paper. The model of an industrial robot built using Autodesk Inventor 3D CAD software was imported into the ANSYS MechanicalANSYS Parametric Design Language (APDL )FEM solver in the "*.STP" format. The material of the robot arm is aluminum 6061 with modulus of elasticity E = 69 GPa, density ρ = 2700 kg/m 3 , and Poisson's ratio ν = 0.3, applied force at the distal-end joint FY = 1300 N, moment MX = 277 Nm, and MZ = 512 Nm.
The 2nd robot arm, together with the subsequent arms and distal-end joint, is shown in Figure  3a. All the subsequent arms are treated as point masses (also known as remote forces) which are used to transfer the force at the distal-end joint (right side) to the 2nd robot arm. Consequently, a resultant force FY and two resultant moments MX and MZ are applied to the 2nd robot arm as shown in Figure  3b. The FEM mesh of the 2nd robot arm is shown in Figure 4a, and the supporting boundary is placed at the proximal end (left side) of the robot arm with the harmonic reducer treated as a fixed connection as shown in Figure 4b. The FEA results are shown in Figure 5.    In this paper we propose a surrogate-based evolutionary optimization (SBEO) method which integrates surrogate models with the multi-objective evolutionary algorithm based on decomposition (MOEA/D-DE). A conceptual sketch is shown in Figure 6. The evolutionary optimization is carried out in the Python programming language using the "PyGMO" open-source optimization library. Herein, we will only give a brief introduction, and the reader can refer to [25,26] for more details.
The procedure of the proposed SBEO method is outlined in the following steps: (1) First, a DOE plan is built based on the Box-Behnken design, and initial sample points are chosen according to the plan matrix. The objective (performance) functions are then analyzed and evaluated using the ANSYS solver following the baseline design CAD model of the robot arm.  4) The weighted sum method, in which each weighting factor is chosen according to the importance or preference of each objective, is then used to determine the best candidate among the generated population samples. These candidates are used to evaluate the fitness with the termination conditions given in Equations (15)- (17). The reduced design space is thereby found to update the surrogate model.  In this paper we propose a surrogate-based evolutionary optimization (SBEO) method which integrates surrogate models with the multi-objective evolutionary algorithm based on decomposition (MOEA/D-DE). A conceptual sketch is shown in Figure 6. The evolutionary optimization is carried out in the Python programming language using the "PyGMO" open-source optimization library. Herein, we will only give a brief introduction, and the reader can refer to [25,26] for more details.
The procedure of the proposed SBEO method is outlined in the following steps: (1) First, a DOE plan is built based on the Box-Behnken design, and initial sample points are chosen according to the plan matrix. The objective (performance) functions are then analyzed and evaluated using the ANSYS solver following the baseline design CAD model of the robot arm.  4) The weighted sum method, in which each weighting factor is chosen according to the importance or preference of each objective, is then used to determine the best candidate among the generated population samples. These candidates are used to evaluate the fitness with the termination conditions given in Equations (15)- (17). The reduced design space is thereby found to update the surrogate model.

Criteria of Model Renewal and Termination
The first termination condition is about the distance between samples and difference of their objective responses: The second one is about the maximum number of affordable expensive function evaluations: In general, if one of the preceding conditions is reached, the optimization is terminated [27].
According to the main effect analysis of DOE in obtaining the sensitivity of design variables, a selected reduced design space is used to renew the RS models. When one of the termination criteria given in Equations (15)-(17) is satisfied, the iterations for model renewing will stop.

Design Procedure of RSM with BBD
The design procedure of the response surface method in association with Box-Behnken design is as follows: (1) Identify the definition and the value of each variable; (2) Create a Box-Behnken design matrix;

Construction of the BBD Matrix
Box-Behnken design is a rotatable second-order design based on three-level incomplete factorial designs, which can decrease the number of experiments to reduce laboratory work. A detailed description of the method can be found in [28]. Herein, a three-factor and three-level Box-Behnken design matrix with 46 runs was chosen as the basic form of the RSM. The Model-3 case is shown in Table 2.

Formulation of the RS Regression Model
The response surface method is used to construct a suitable approximation of the true functional relationship between a dependent variable (the response) and a vector of independent variables. The response is generally evaluated experimentally or analyzed using numerical simulation [29].
A quadratic regressive equation was chosen to be the model of the robot arm, as shown below: where Y is the response or performance function, x i , x j are the decision variables, β o is the offset term, β i represents the linear term, β ii denotes the quadratic term, β ij denotes the interactive term, and ε is the statistical error. The relationship between the coded and real variables is described in Equation (8). The values of β are obtained from Equation (19) using inverse matrix transformation by Python's numpy library's "Polynomial.fit" function as below:

Calculation of the Coefficients of the RS Model
The coefficients of different surrogate models obtained via Equation (18) and the real variables obtained by the transformation of the coded variables are shown in Table 3. Table 3. The coefficients of RS models and real variables.

Validation of the RS Model
In statistical analysis, the coefficient determination of R 2 [30] is a measure used to assess how well a model can explain and predict outcomes. It is also commonly known as "R-squared", to be used as a guideline for measuring the accuracy of the model. It ranges from 0 to 1 and will be closer to 1 when the model is more precise. Joglekar and May [31] suggested that the value of R 2 for the best-fitting model should be at least 0.8. According to each DOE matrix, the obtained model is able to produce a predicted response (called Y_Model). The analyzed response produced by the CAE solver is called Y Exact. If these two values are equal, the red points in Figure 7 will lie on the diagonal line. The R 2 values of W, J, and D models being, respectively, 0.9875, 0.9865, and 0.9229 mean that the surrogate models evaluated in this study can be used to explain the W, J, and D very well.

Validation of the RS Model
In statistical analysis, the coefficient determination of R 2 [30] is a measure used to assess how well a model can explain and predict outcomes. It is also commonly known as "R-squared", to be used as a guideline for measuring the accuracy of the model. It ranges from 0 to 1 and will be closer to 1 when the model is more precise. Joglekar and May [31] suggested that the value of R 2 for the bestfitting model should be at least 0.8. According to each DOE matrix, the obtained model is able to produce a predicted response (called Y_Model). The analyzed response produced by the CAE solver is called Y Exact. If these two values are equal, the red points in Figure 7 will lie on the diagonal line. The R 2 values of W, J, and D models being, respectively, 0.9875, 0.9865, and 0.9229 mean that the surrogate models evaluated in this study can be used to explain the W, J, and D very well.

Multi-objective Evolutionary Optimization
The use of MOEA/D-DE algorithm for solving multi-objective optimization in robot arm design allows us to find the Pareto front and the corresponding optimal parameter sets. The surrogate models in Table 2 were substituted into the MOEA, MOEA/D-DE operators with the adoption of the following parameters: (1) Generations = 300; (2) Weight generation = grid → method to generate the weights; (3) Decomposition method = Tchebycheff → method used to decompose the objectives; (4) Neighbor = 20 → size of weights neighborhood; (5) Crossover parameter = 0.9; (6) Differential evolution parameter = 0.5; (7) Distribution index = 20 → used for polynomial mutation; (8) Real b = 0.9 → chance that the neighborhood is considered at each generation, rather than the whole population (part of the diversity preservation mechanism); (9) Limit = 2 → number of copies reinserted in the population (part of the diversity preservation mechanism); (10) Preserve diversity = True → diversity preservation mechanisms.
We also set the preference criteria shown in Equation (5)

Multi-objective Evolutionary Optimization
The use of MOEA/D-DE algorithm for solving multi-objective optimization in robot arm design allows us to find the Pareto front and the corresponding optimal parameter sets. The surrogate models in Table 2 were substituted into the MOEA, MOEA/D-DE operators with the adoption of the following parameters: 1) Generations = 300; 2) Weight generation = grid  method to generate the weights; 3) Decomposition method = Tchebycheff  method used to decompose the objectives; 4) Neighbor = 20  size of weights neighborhood; 5) Crossover parameter = 0.9; 6) Differential evolution parameter = 0.5; 7) Distribution index = 20  used for polynomial mutation; 8) Real b = 0.9  chance that the neighborhood is considered at each generation, rather than the whole population (part of the diversity preservation mechanism); 9) Limit = 2  number of copies reinserted in the population (part of the diversity preservation mechanism); 10) Preserve diversity = True  diversity preservation mechanisms.
We also set the preference criteria shown in Equation (5)

Simulation Results
The optimization problem stated in Equations (1)-(4) for the optimal shape design of the 2nd robot arm was solved using the proposed method. The results related to the optimal design are given in the following to illustrate the effectiveness and feasibility of the proposed method.

Simulation Results
The optimization problem stated in Equations (1)-(4) for the optimal shape design of the 2nd robot arm was solved using the proposed method. The results related to the optimal design are given in the following to illustrate the effectiveness and feasibility of the proposed method.

Pareto Front History
In solving a multi-objective optimization problem, a solution that has one or several objective functions that are impossible to further optimize while other objective functions are not deteriorating is called "Pareto optimal" [32]. In this case, the Pareto front history of MOEA is plotted in Figure 9, and it is obvious that the solution started to converge at the 10th generation of MOEA.

Pareto Front History
In solving a multi-objective optimization problem, a solution that has one or several objective functions that are impossible to further optimize while other objective functions are not deteriorating is called "Pareto optimal" [32]. In this case, the Pareto front history of MOEA is plotted in Figure 9, and it is obvious that the solution started to converge at the 10th generation of MOEA.

Weighted Sum History
The multi-objective optimization problem of Equation (5) was evaluated to obtain the weighted sum result. The objectives and design variables in different models and generations according to the preference criteria on the optimal candidate in each generation of MOEA are shown in Figure 10a and

Weighted Sum History
The multi-objective optimization problem of Equation (5) was evaluated to obtain the weighted sum result. The objectives and design variables in different models and generations according to the preference criteria on the optimal candidate in each generation of MOEA are shown in Figures 10a  and 10b. Optimal values for W of 33.02 kg, J of 6.3328 kg.m 2 , and D of 0.0491 mm were obtained for H i = 198.81 mm, H m = 549.99 mm, V m = 100.07 mm, R h = 60.00 mm, and T a = 90.13 mm.

Pareto Front History
In solving a multi-objective optimization problem, a solution that has one or several objective functions that are impossible to further optimize while other objective functions are not deteriorating is called "Pareto optimal" [32]. In this case, the Pareto front history of MOEA is plotted in Figure 9, and it is obvious that the solution started to converge at the 10th generation of MOEA.

Weighted Sum History
The multi-objective optimization problem of Equation (5) was evaluated to obtain the weighted sum result. The objectives and design variables in different models and generations according to the preference criteria on the optimal candidate in each generation of MOEA are shown in Figure 10a and

The Improvement of Performance
The predicted optimal parameter sets for the models of the 2nd robot arm obtained in the SBEO process were analyzed using ANSYS to obtain the responses of the selected objectives, which were then compared with those of the baseline design. The results of different models are shown in Table 4 for comparison. According to the data, the first design (Model-1) greatly improved the deformation objective, but showed no significant improvements for the other objectives. In the renewal design Model-3, the objective W was decreased by 16.1%, J was decreased by 22.8%, and D was decreased by 20.3%; optimal values of W = 33.08 kg, J = 6.4906 kg.m 2 , and D = 0.0513 mm were obtained for H i = 198.81 mm, H m = 549.99 mm, V m = 100.07 mm, R h = 60.00 mm, and T a = 90.13 mm. The results reveal that our proposed method is able to solve this real engineering optimization problem, and its effectiveness is thus validated.

The Improvement of Computational Efficiency
A computational efficiency comparison between the use of surrogate models and the use of ANSYS is shown in Table 5. If ANSYS is used, the computational time for a single configuration is 30 s. We set the computation of each configuration as a sequential process (computing one configuration at a time) performed in a single computer (no cluster computing). It is feasible to compare the computational cost of each model. The results of the comparison show that the computational cost of the proposed method is significantly less than that of ANSYS.

The Relationship between Objectives and Design Variables
There are two ways to represent the relationship between design variables and objective functions: main effect analysis based on the results of design of experiment and pairwise plots of the SBEO results.

Main Effect Analysis
A main effect is frequently used in the context of factorial designs and regression models to distinguish main effects from interaction effects. A main effect occurs when the mean response changes across the levels of a factor. According to the ranges of design variables in Table 1, the lower bound is denoted as level (−1), the upper bound is denoted as level (+1), and the average of the two is denoted as level (0) when running the DOE. The main effect plots for the objectives of weight (W), deformation (D), and moment of inertia (J) are shown in Figure 11. Several results observed in the main effect plots are stated below.
On the objective W: (1) All factors, except R h , can cause the response to increase when they move from low level (−1) to high level (+1). While H i and H m can only cause the response to increase slightly as the level changes, V m and T a can cause the response to increase significantly, especially T a . (2) Factor R h can cause the response to decrease when its level increases, but the slope of the response only changes slightly. (3) Consequently, all five factors appear to have a main effect on the response W. Among them, T a is the most important factor because of its high significance.
On the objective J: (1) Factors H i and R h present no significant main effects because the response remains nearly the same when they move from low level to high level. (2) Factors V m and T a can cause the response to increase significantly when they move from low level to high level, and the slope of the increased response induced by T a is the sharpest. On the objective D: (1) When the factors move from low level to high level, only R h causes response D to decrease, while the others cause response D to increase. (2) Factors H i , V m , and T a cause the response D to decrease sharply when the factors move from low level to high level, and the slopes of the responses remain almost the same. (1) All factors, except Rh, can cause the response to increase when they move from low level (−1) to high level (+1). While Hi and Hm can only cause the response to increase slightly as the level changes, Vm and Ta can cause the response to increase significantly, especially Ta. (2) Factor Rh can cause the response to decrease when its level increases, but the slope of the response only changes slightly. (3) Consequently, all five factors appear to have a main effect on the response W. Among them, Ta is the most important factor because of its high significance.
On the objective J: (1) Factors Hi and Rh present no significant main effects because the response remains nearly the same when they move from low level to high level. (2) Factors Vm and Ta can cause the response to increase significantly when they move from low level to high level, and the slope of the increased response induced by Ta is the sharpest. On the objective D: (1) When the factors move from low level to high level, only Rh causes response D to decrease, while the others cause response D to increase. (2) Factors Hi, Vm, and Ta cause the response D to decrease sharply when the factors move from low level to high level, and the slopes of the responses remain almost the same.

Pairwise Plot Analysis
A pairwise plot of the relations between the objectives and design variables was obtained and is shown in Figure 12. At the first generation, we can see that our evolutionary algorithm takes random samples for each variable from x1 to x5 within the defined lower and upper bounds. Using Kernel Density Estimation (KDE), we can visualize the distribution for each parameter with respect to other variables and confirm that the sampling was indeed randomly distributed. Histogram plots also confirm that the initial population was distributed throughout the parameter space (lower and upper

Pairwise Plot Analysis
A pairwise plot of the relations between the objectives and design variables was obtained and is shown in Figure 12. At the first generation, we can see that our evolutionary algorithm takes random samples for each variable from x 1 to x 5 within the defined lower and upper bounds. Using Kernel Density Estimation (KDE), we can visualize the distribution for each parameter with respect to other variables and confirm that the sampling was indeed randomly distributed. Histogram plots also confirm that the initial population was distributed throughout the parameter space (lower and upper bounds of each variable). Interestingly, the scatter plots and correlation coefficients reveal that the proposed mathematical model can be attained in the parameter space of variables. Only variable x 5 is linearly correlated with the objectives W, D, and J. Also, W, D, and J are linearly correlated with each other. In detail, W has a negative correlation with D, and D has a negative correlation with J. However, W and J have a positive correlation. So, we can conclude that for smaller values of W and J, we can expect to have a larger value of D. Similarly, for a smaller value of D, we expect to have larger values of W and J. From a mechanics point of view, this can be seen as (1) W being directly proportional to J, (2) W being inversely proportional to D, and (3) J being inversely proportional to D. Thus, we can conclude that our surrogate model captures the mathematical properties of the system nicely.

Discussion
In this work, a surrogate-based evolutionary optimization method together with Pareto front was presented to design the 2nd robot arm of an industrial robot for attaining high speed performance of the robot. The robot arm was optimized considering the minimum weight, deformation, and moment of inertia. The optimization problem was subjected to constraints on performance indicators such as the W2SR, W2JR, and geometric conditions. Five design variables selected from the important shape geometric parameters were placed under realistic engineering consideration in the optimal design. In order to solve the multi-objective problem, the weighted sum method was used to determine the best candidate from the generated population samples, where the weighting factors were chosen according to the importance or preference of each objective. The characteristics of the optimized 2nd robot arm were also compared with those of the baseline 2nd robot arm. The results obtained in the comparison show that the optimized profile not only increased the performance indicators W, J, and D by up to 16.1%, 22.8%, and 20.3%, respectively, but also saved over 90% on the computational cost. These advantages indicate the effectiveness of the proposed SBEO method when compared with the existing ANSYS-based methods and thus confirm the efficiency of the proposed method in solving complex realistic engineering problems. Some conclusions can be drawn as follows:

Discussion
In this work, a surrogate-based evolutionary optimization method together with Pareto front was presented to design the 2nd robot arm of an industrial robot for attaining high speed performance of the robot. The robot arm was optimized considering the minimum weight, deformation, and moment of inertia. The optimization problem was subjected to constraints on performance indicators such as the W2SR, W2JR, and geometric conditions. Five design variables selected from the important shape geometric parameters were placed under realistic engineering consideration in the optimal design. In order to solve the multi-objective problem, the weighted sum method was used to determine the best candidate from the generated population samples, where the weighting factors were chosen according to the importance or preference of each objective. The characteristics of the optimized 2nd robot arm were also compared with those of the baseline 2nd robot arm. The results obtained in the comparison show that the optimized profile not only increased the performance indicators W, J, and D by up to 16.1%, 22.8%, and 20.3%, respectively, but also saved over 90% on the computational cost. These advantages indicate the effectiveness of the proposed SBEO method when compared with the existing ANSYS-based methods and thus confirm the efficiency of the proposed method in solving complex realistic engineering problems. Some conclusions can be drawn as follows: (1) The optimal shape design of a robot arm performed in this study showed that shape design optimization is useful to improve the performance of industrial robots with real engineering considerations. (2) The optimal design of the 2nd robot arm obtained via the multi-objective optimization approach features good performance for conflicting objectives, and a Pareto front was obtained by employing preference sets of weighted factors for the various design objectives. (3) Compared with the conventional expensive ANSYS-based simulation technique, the proposed SBEO is more practicable and effective for shape optimization of robot arms parameterized with many design variables. In the SBEO, the population samples and Pareto fronts are able to converge quickly to the optimized hot zone of design candidates according to the designer's preference criteria to obtain better performance. (4) It was shown that objective W is proportional to J, and both W and J are inversely proportional to D. The objectives are all affected by the five adopted design variables, of which thickness T a has the most significant effects, which depend on the selected criteria. When T a , R h , and V m are fixed, smaller H i and H m at the same time could cause larger deformation in the robot arm design optimization.
For future research directions, it is worthwhile to seek further applications of the SBEO algorithm to robot design, as well as other types of robot mechanism. In this study, we only considered the shape design of the 2nd robot arm, but did not consider the design of the whole body and robot dynamic behavior yet. Anyone who is interested in this area can extend our study to the abovementioned topics based on our results. In addition, an experimental study to test and validate the performance of the algorithm is also recommended.