Dimensional Variation Analysis for Rigid Part Assembly With an Improvement of Monte Carlo Simulation

Dimensional variation has significant effect on the quality of product. Recently, Monte Carlo (MC) simulation is widely used in dimensional variation analysis, with high accuracy and adaptability, but there is the problem of low computational efficiency. Aiming to address this problem, an improvement of MC simulation is proposed through a two-phases solution. In the first phase, surrogate model is used to approximate the locating constraint equations for a 3D part, which reduces the nonlinear coupling between dimensional deviations and avoids large-scale solution of nonlinear equations in dimensional variation analysis on the condition of ensuring the accuracy. In the second phase, random samples used in MC simulation are replaced by low discrepancy sequences, which enable the samples to have better homogeneity and representativeness and reduce the number of samples required in the dimensional variation analysis. Finally, two examples are used to demonstrate the effectiveness of the method, and the results show that the two-phases solution has advantages both in the accuracy and efficiency.


I. INTRODUCTION
Dimensional variation is a significant index for evaluating product's quality. Oversize dimensional variations may influence the product appearance and operational performance, and it may even result in excessive assembly stress, shorten the fatigue life of products, and impact the security of products. Variations may be unavoidable in manufacturing, and it may cause variations in geometrical features. Meanwhile, these variations may transfer and accumulate gradually with the assembly processes, and as a result, the assembly dimensions may exceed the specification. According to the statistics of pertinent literature, two thirds of product's quality problems are caused by the dimensional variations [1]. For the reason, engineers need to know the influence of dimensional variations on the quality of product. Thus, the method of dimensional variation analysis is necessary to answer the question. Dimensional variation propagation model (DVPM) is fundamental for dimensional variation analysis, which The associate editor coordinating the review of this manuscript and approving it for publication was Ming Luo . is used to describe the accumulation of deviations in the assembly. Several DVPMs are developed during past decades, including matrix model [2], [3], Jacobian-Torsor model [4], [5], T-Map model [6], [7], DLM model [8], [9], stream-of-variation (SOV) model [10]- [12] and skin shapes mode [13]- [15]. A summary of the literature is presented as follows: • The matrix model is developed from the concept of Topological and Technological Related Surface (TTRS), in which, dimensional deviations and their accumulations are represented by a series of homogeneous transformation matrixes. The theoretical basis of matrix model is linear operation of homogeneous transformation matrixes.
• In the Jacobian-Torsor model, dimensional deviations and their accumulations are represented by small displacement torsor (SDT) and Jacobian, respectively. Due to the complexity of parallel mechanism or structure in product, it is hard to analyze the accumulations of dimensional deviations in parallel mechanism or structure, which may result in a non-linear accumulation of dimensional deviations.
• The T-Map model uses convex set to represent the interval of dimensional deviations, and calculates the accumulation of dimensional deviations through Minkowski sum of the convex set. A specific program is necessary for T-Map model due to the complex mathematical representation of convex set and Minkowski sum.
• The DLM model considers dimensional deviations as independent vectors, and the accumulation of dimensional deviations is derived from the assembly constraints, which is represented by an interconnection of the vectors. For complex product structure, it is not easy to generate a DLM model.
• The SOV model is based on the theory of linear state space, in which, dimensional deviations and their accumulations are represented by state vectors and state transfer functions, respectively. It is worth pointing out that the dimensional deviations in part-to-part joints and part-to-fixture joints are both considered in SOV model, which is suitable for complex assemblies.
• The skin shapes model is recently developed DVPM that reflects the geometric deviations by discretizing geometric features into points and imparting different geometric deviations to the points. The difficulty is to ensure that the geometric deviations of different discrete points are within tolerances. Recently, the skin shapes model has been applied in the field of tolerance analysis with consideration of many kinds of mechanical joints. With a summarization of the DVPMs in the literature, we can find that there exist many models or methods in industry and academia. However, common shortcoming needs to be deeply analyzed. For most of the DVPMs, it is mainly assumed that the dimensional deviations have linear relationship or that the relationship between dimensional deviations is simplified as linearity. In fact, due to the complexity of products and fixtures, parts in an assembly process will be constrained by multiple mating constraints, which may result in multiple paths and nonlinearities in the accumulation of dimensional deviations. That is to say, the relationship between dimensional deviations is non-linear. In this paper, we will explain this non-linear relationship in detail.
The purpose of dimensional variation analysis is to determine the impact of the accumulation of dimensional deviations on product quality based on DVPM. In general, the outputs of dimensional variation analysis are statistical characteristics of a functional dimension (i.e. a gap or an assembly dimension) in product, such as mean value and standard deviation of the functional dimension. Currently, typical methods for dimensional variation analysis include worst-case method, root sum square (RSS) method, statistical analysis methods and Monte Carlo method [16]. A summary of the literature is presented as follows: • The worst-case method is used to analyze the conditions when the value of dimension reach maximum or minimum [17], [18]. It is rarely used in engineering because results are usually too much pessimistic.
• In the RSS method, the dimensional deviations in assembly are assumed to comply with certain statistical distributions, and the dimensional variations in assembly are calculated based statistics theory. According to literature, the results of RSS method are usually too much optimistic [19], [20].
• The Monte Carlo method is the most widely used method, especially applicable for complex assembly. Several commercial tolerance analysis software such as CETOL, eM-TolMate, VisVSA and 3DCS all use Monte Carlo method [21].  [32], and Liu et al. [33]. With analysis of the research, the Monte Carlo method is the most widely used in dimensional variation analysis, with high accuracy and adaptability. It is found that the theoretical foundations of Monte Carlo method are the law of large numbers and the central limit theorem. These two laws theoretically ensure that as long as enough samples of dimensional deviations are generated, the result of assembly variation analysis will possess high accuracy and reliability. Thus, the accuracy of Monte Carlo method is related to the number of samples. In actual application, lots of samples are required to improve the accuracy, which spends too much time. It is pointed that in order to improve the accuracy of MC method with 10 times, the number of samples is needed to be increased for nearly 100 times [31]. For the reason, the efficiency of Monte Carlo method is an important factor that affects its further application in engineering, which is exactly the purpose of this study.
The paper is organized as follows. In section II, the preliminary background and comprehensive conceptual solution is represented, in which, the problem formulation and the overview of the proposed the two-phases solution are introduced. The proposed method, including two phases, is presented in section III. In section IV, two cases are studied, including an illustrative example and an industrial case. Finally, conclusions are represented in section V.

II. PRELIMINARY BACKGROUND AND COMPREHENSIVE CONCEPTUAL SOLUTION A. DIMENSIONAL VARIATION PROPAGATION MODELING
Without loss of generality, dimensional variation propagation process for an assembly station can be expressed with Figure 1. In this paper, we only consider single station assemblies, which do not consider the dimensional variation propagation between multi-stations. However, it is noted that the proposed method can be extend to multi-stations.  As illustrated in Figure 1, part k is influenced by three kinds of deviations, which are presented by h k−1 , e k and ξ k , respectively. h k−1 represents the deviation transferred from part k-1; e k represents the deviation transferred from the fixture of part k; ξ k represents random deviation. Observation point can be selected on part k to observe or measure dimensional variations in the assembly, which is presented by δ k in Figure 1.
In actual assembly, parts are located by assembly constraints, such as against, fit, etc [26]. That is to say, assembly constraints are the carrier and bridge for dimensional variation propagation. According to Huang et al. [10], [11] and Tang et al. [34], assembly constraints applied to a part can be transformed into six points, which are referred to as locating points (LPs). The six LPs can be considered as the carrier of variations in parts and fixtures, which are used to establish the relationship between deviations in the assembly. In actual application, the following principles are adopted to select the six LPs. Firstly, three LPs are selected on the primary datum surface, which make the area as large as possible. Secondly, two LPs are selected on the secondary datum surface, which make the distance as long as possible. Thirdly, a LP is selected in center of the tertiary datum surface. As illustrated in Figure 2, the block is located by six LPs. L 1 , L 2 and L 3 are selected in the primary datum surface. L 4 and L 5 are selected in the secondary datum surface. L 6 is selected in the tertiary datum surface. With these six LPs, the block is fully constrained in the assembly.
It is pointed out that the foregoing simplification of the six LPs holds only for rigid parts assembly. Based on the simplification, the assembly of a part can be transformed into a general locating layout shown in Figure 3, which is used to establish the relationship between deviations in assembly. In Figure 3, the six LPs are represented as L i , where i = 1, 2, · · · , 6. OXYZ is the reference coordinate system, and O X Y Z is the local coordinate system attached to the part; r i is the positional vector of L i in OXYZ ; r i is the positional vector of L i in O X Y Z . Below, unless otherwise specified, the vector with is described in O X Y Z . When there are deviations in LPs, the spatial position of the part will change, and the core issue of assembly variation analysis is to determine the impact of LPs' deviations on the spatial position of the part.
In actual assembly, the six LPs should contact with the part to realize the function of location. Based on this actual condition, it is assumed that the surface of the part in contact with L i is S i (x , y , z ) = 0, and L i is always in contact with S i (x , y , z ) = 0. For complex product, S i (x , y , z ) = 0 is probably a freeform surface. The second-order Taylor expansion of S i (x , y , z ) = 0 at r i in Figure 3 can be expressed as follows.
n ixx n ixy n ixz n iyx n iyy n iyz n izx n izy n izz   (2) By translating r into the reference coordinate system (OXYZ ), the expression (1) can be represented as follows.
In which, i (ψ, r i ) represents the locating constraint equation of L i ; ψ = r T 0 , T T represents the six posture parameters; r 0 = [r 0x , r 0y , r 0z ] T represents the positional vector of the origin point of O X Y Z in OXYZ ; where c and s are the abbreviation of cos and sin, respectively. Without loss of generality, we assume O X Y Z coincides with OXYZ when the part is in its original position. That is to say, when the part is in the original position, there exists r i = r i and ψ = 0. Based on the above assumption, when there are deviations in L i , we can derive the vector of r i as follows.
In which, δ x , δ y and δ z represent the deviations of observation point in the axes of x, y and z, respectively; δ x , δ y and δ z respectively represent the coordination of observation point in the axes of x, y and z. It can be learnt that the solving of expression (6) is the key issue of assembly variation analysis. However, it is not difficult to see that expression (6) is an equation set consisting of six nonlinear equations, and each equation is implicit. Therefore, it is difficult to find the solution of expression (6). On the one hand, numerical method will be used to solve the expression (6), which is very time-consuming. On the other hand, in assembly variation analysis, engineers also need to obtain the statistical characteristics of observation point, such as the mean and the standard deviation. However, the statistical characteristics of observation point may not be calculated in an analytical way due to the high nonlinearity of expression (6). Therefore, engineers usually adopt Monte Carlo (MC) method to obtain statistical characteristics of the observation point. It will cost engineers' lots of time, especially when the number of parts in a product is too large. The difficult point noted above is exactly the starting point of this paper. In the following, we will propose a two-phases solution to help engineers complete assembly variation analysis with high precision and high efficiency.

B. OVERVIEW OF THE PROPOSED TWO-PHASES SOLUTION
The two-phases solution proposed in this paper includes two parts: surrogate modeling and Quasi-Monte Carlo (QMC) simulation. The general idea of the two-phases solution is shown in Figure 4.
Phase I: Firstly, several optimized samples of LPs' deviations are generated. Secondly, a numerical method is used to calculate the posture parameters corresponding to each sample, and we define these posture parameters as responses. Thirdly, by combining with the samples of LPs' deviations and the corresponding responses, a high precision approximation relationship of expression (6) is constructed. The approximation relationship transforms the relationship between LPs' deviations and the responses from implicit into explicit, and it reduces the nonlinearity of expression (6) at the same time.
Phase II: Aiming at the limitation of MC method that it needs a large number of random samples, QMC simulation is proposed. In the QMC simulation, low discrepancy sequences (LDS) is adopted to replace random samples, which can effectively find the statistical characteristics of the observation point.
By combining Phase I and Phase II, we can avoid largescale solving of nonlinear equations and greatly improve the efficiency of assembly variation analysis.
According to the general idea of the two-phases solution presented in Figure 4, the procedure of the two-phases solution is proposed in Figure 5.

III. THE PROPOSED METHOD
A. PHASE I: SURROGATE MODELING Surrogate modeling plays a very important role in two-phases solution and determines the accuracy of assembly variation VOLUME 8, 2020 analysis. The overall procedures of surrogate modeling include two steps [35]. Firstly, obtain a certain number of exact inputs and outputs (or responses) of expression (6). Then, choose a reasonable approximation model to fit the relationship between exact inputs and outputs, which has low complexity and high efficiency. Kriging method is a commonly used method in surrogate modeling, which has fine fitting accuracy for non-linear relationship [36]. We will choose Kriging method for surrogate modeling.
In order to get exact inputs and outputs of expression (6), design of experiment (DOE) is needed to get optimized samples of LPs' deviations. Latin hypercube sampling (LHS) is selected for DOE, which has high efficiency and good uniformity [35], [37]. In general, it is assumed that the deviations of L i obey normal distribution, and the mean and standard deviation of the normal distribution are represented by µ i and σ i , respectively. We use the functions provided by MAT-LAB R to generate LHS according the normal distribution. In expression (5), the deviations of L i is represented by w i , so we define the LHS of w i as 1 ω i , 2 ω i , · · · , j ω i , · · · , m ω i , in which, j represents the jth LHS; m represents the number of LHS; j = 1, 2, · · · , m. For each LHS, the equation (6) can be rewritten as: where, j r i represents the vector of L i when the deviation of L i is j ω i . It is defined that j r = [ j r 1 T , j r 2 T , · · · , j r 6 T ] T and j ψ. Where, j ψ represents the posture parameters when the vectors of LPs are j r. Substitute j r and j ψ into equation (6), and then we can obtain the following expression. j ψ, j r = 0 Expression (9) is in fact an equation set consisting of six equations, which is the instantiation of expression (6). Levenberg-Marquardt (LM) is adopted to find the solution of j ψ, and the detailed procedures of LM can be seen in the reference.
According to expression (8) and expression (9), we can know that the six posture parameters are actually determined by LPs' deviations, which are represented by w i in expression (5). Here, it is assumed that the functional relationship between the six posture parameters and LPs' deviations can be expressed as: In which, w = [w 1 , w 2 , · · · , w 6 ] T . The aim of surrogate modeling is to generate an approximation relationship of expression (10) with usage of j ω i (i = 1, 2, · · · , 6) and j ψ. Based on Kriging method, expression (10) can be expressed as follows.
In which, f (w) is a known regression model, usually choosing linear or quadratic polynomial function; γ is the unknown regression coefficient; f T (w) γ constitutes a deterministic part to express the global characteristics of ψ (w); P (w) is a statistical process to express the local characteristics of ψ (w).
In this paper, DACE toolbox of MATLAB R is used to find γ and P (w).

B. PHASE II: QUASI-MONTE CARLO SIMULATION
QMC simulation is used to obtain the statistical characteristics of the observation point. The main difference between QMC and MC is that low discrepancy sequences (LDS) is used to replace the random sample in QMC. For two dimensional variables with range of 0 to 1, one hundred random samples and LDS are illustrated in Figure 6 and 7, respectively. As shown in Figure 6, the distribution of random samples is disordered. In some areas, the samples are too concentrated and even overlap, while in other areas the number of samples is very small and even there is no sample. As shown in Figure 7, the distribution of LDS is very uniform, and there is no overlapping or blank area. Thus, by analyzing the samples from the perspective of their uniformity, the following conclusions can be drawn. Firstly, redundancy is easily to appear in random samples, i.e. in certain area of the sample space, there may appear two or more samples, which is an important reason of low efficiency of MC method. At the same time, due to the randomness of the random samples, the samples are distributed randomly in the sampling space, which will lead to the ''oscillation'' phenomenon in the results of MC method. The ''oscillation'' phenomenon will be illustrated in the case study. Secondly, LDS are uniformly distributed in the sampling space. For example, there will not appear two or more samples in each area of the sample space.
In other words, LDS can be effectively used, which will help improve the efficiency of QMC simulation. There are several methods to generate LDS, i.e. Hua-Wang (H-W) points and Halton sequences. The methods for generating LDS are compared in the reference [38]. In this paper, H-W points will be used as LDS, which can be generated by expression (12).
· · · , mod 2 j cos 2π l p , 1 ] (12) In which, l is the dimension of variables, and we set l = 6 in this paper; p is prime number and p ≥ 2l + 3; n is number of samples; mod represents modular operation. It needs to be note that the samples generated by the calculation of expression (12) are in uniform distribution, which can be transformed to other distributions through a transformation [39].

IV. CASE STUDIES A. CASE 1: AN ILLUSTRATIVE EXAMPLE FOR THE COMPARISON OF DIFFERENT METHODS
As is shown in Figure 8, an ellipsoidal part is selected as example to verify the effectiveness of the proposed method. The ellipsoid is represented by  In order to verify the effectiveness of the proposed method, we will adopt MC method and the two-phases solution to calculate the posture parameters of the ellipsoid. The results of MC method will be taken as the exact solution, and the relative error compared with MC method is expressed as: In which, η represents the results of two-phases method, and η MC represents the result of MC method. The programs for this case are formulated in MATLAB, and the computational environment includes Windows XP, CPU with 2.27HZ and RAM with 1.92GB.

1) THE VALIDATION OF ACCURACY
It is assumed that the deviations of the six LPs (L 1 -L 6 ) are all in normal distributions and the standard deviations of the normal distributions are equal. We adopt the MC method and two-phases solution to calculate the posture parameters, respectively. In the process of calculation, the number of samples used in MC method and QMC simulation are both 25,000, and 60 samples are used in surrogate modeling. We choose two posture parameters (φ and r 0x ) as the object for analyzing, and the results are shown in Figure 9 to 12. As shown in Figure 9, the results of the mean value, calculated by the MC method and the two-phases solution, do not fit very well. The relative errors are large. However, as we can see that the mean values are mainly between 10 −5 and 10 4 in Figure 9. In actual engineering, the tolerances are difficult to have such high requirements. Meanwhile, the mean values in Figure 9 are infinitesimal when they are comparing to Figure 10, 11 and 12. Thus, the relative errors in Figure 9 may be caused by some random perturbations when the methods are sampling. We believe that the MC method and the two-phases solution both have high accuracy when they are calculating the results in Figure 9.
As shown in the Figure 10 to 12, the results of two-phases solution are almost the same as MC method, and the relative errors are no more than 1%. VOLUME 8, 2020   Therefore, we can draw the conclusion that the results of two-phases solution are almost consistent with MC method, which has high accuracy.

2) THE VALIDATION OF EFFICIENCY
The two-phases solution and MC method have similar accuracy. In order to compare the efficiency of these two methods, the relationship between the standard deviation of posture parameter φ and the number of sample is analyzed. It is   assumed that the standard deviations of LPs are 0.5, and then we get the relationship between the standard deviation of posture parameter φ and the number of sample, as illustrated in Figure 13. As can be seen from the figure, two-phases solution and MC method both show instability when the number of samples is small. However, the convergence rate of two-phases solution is greatly faster than MC method along with the increasing of the number of samples, especially when 5868 VOLUME 8, 2020 the number of sample reaches 6000. At the same time, the convergence of two-phases solution is stable, while the convergence of MC method has an ''oscillation'' phenomenon. Thus, two-phases solution not only inherits the accuracy of MC method, but also improves the convergence rate and stability. Figure 14 is the relationship between the number of samples and computing time corresponding to Figure 13. In overall, the computing time of the two-phases solution is much shorter than MC method. For example, when the number of samples is 5000, the computing time of the two-phases solution is 1.109 seconds, while the computing time of MC method is 50.187 seconds; when the number of samples is 25000, the computing time of the two-phases solution increases to 3.297 seconds, while the computing time of MC method increases to 255.190 seconds. The time of MC method is nearly 100 times bigger than the time of two-phases solution. The reason why two-phases solution has a higher efficiency is that the application of surrogate model greatly reduces the time of solving nonlinear equations.

B. CASE 2: A SIMPLIFIED FUSELAGE
A simplified fuselage is shown in Figure 15a. The fuselage consists of three parts including left frame (X-00-02-01), right frame (X-00-02-02) and skin (X-00-02-03). The selection of the LPs on the right frame and skin is shown in Figure 15b and Figure 15c. In the figure, the direction of the arrow represents the normal vector of LP, and the endpoint of the arrow indicates the positional vector of LP. The distribution of LPs on the left frame is similar with that of the right frame, so its schematic figure is not shown. As is shown in Figure 15d, six observation points (JS1 -JS6) on the skin are selected to observe the contour of the fuselage. The positional vectors and normal vectors of the LPs and the observation points are seen in the appendix.

1) VARIATION ANALYSIS
As illustrated in Figure 16, the dimensional variation propagation process of the fuselage includes three steps. Firstly, the left frame is assembled, which is located by the six points, X-00-02-01-L1 to X-00-02-01-L6. With analysis of  the general locating layout of the left frame, the dimensional variations of the locating point, X-00-02-03-L1, can be obtained. It should be pointed that X-00-02-03-L1 is not only a point on the left frame, but also a point on the skin. That is to say, X-00-02-01-L1 is derived from the contact constraint between the left frame and the skin, which bridges the left frame and the skin. Secondly, with analysis of the general locating layout of the right frame, the dimensional variations of the locating points, X-00-02-03-L2 and X-00-02-03-L3, can be obtained, respectively. Thirdly, the skin is located by the six points, X-00-02-03-L1 to X-00-02-03-L6. With analysis of the general locating layout of the skin, the dimensional variations of the observation points can be calculated.
According to the process data of an aircraft manufacturing enterprise, it is assumed that the standard deviations of LPs on the left and right frames are in normal distributions. The deviations of LPs on the left obey a normal distribution N (0, 0.0883), and the deviations of LPs on the right frames obey a normal distribution N (0, 0.1667). In the stage of product design, the tolerances of the six observation points are set to ±0.6mm. The MC method and two-phases solution are adopted to analyze the standard deviations of the observation points. The results are presented as Table 1.

2) ANALYSIS OF THE RESULTS
The results in Table 1 are transformed into a bar chart to illustrate the difference of the two methods, as shown in Figure 17. It can be seen from the figure that the results of MC method and two-phases solution are very close, and their relative error is no more than 2.5%. Therefore, the proposed two-phases solution meet the requirements of industrial usage.
In the stage of product design, engineers need to know whether the tolerance of observation points is in design value  or not. The tolerance, which is nominal in design stage, is ±0.6mm. According the principle of 6-sigma, the tolerance of each observation point is equal to 6sigma. Here, Sigma is the standard deviation of observation point. That is to say, the standard deviations of the six observation points need to less than 0.2mm, Otherwise, it is out of design. As illustrated in Figure 18, both the MC method and two-phases solution indicate that the standard deviation of JS4 is out of design. Thus, a procedure named tolerance optimization is necessary to find feasible deviations of LPs.

V. CONCLUSION
Dimensional variation analysis for assembly plays an important role in product design and manufacturing, which can help engineers carry out tolerance analysis, tolerance optimization and fault diagnosis. This paper proposes a two-phases solution of dimensional variation analysis. In the first phase, the approximate relationship between dimensional deviations is established based on surrogate modeling, which makes the relationship between dimensional deviations transform from implicit into explicit and reduces the nonlinear coupling between dimensional deviations. In the second phase, QMC is proposed, in which, LDS replace random samples to reduce the number of samples in variation analysis. Combining respective advantages of the two phases, the problem of low efficiency in MC method is improved, and the advantage of MC method, high accuracy, is also inherited.
In the subsequent research, we will primarily study the following questions. Firstly, efficient software for dimensional variation analysis will be studied, which will integrate the proposed method into current CAD system. Secondly, methods for tolerance analysis, tolerance optimization and fault diagnosis will be studied.

APPENDIX
The positional vectors and normal vectors of the LPs and the observation points in case 2 are presented in Table 2.
WENBIN TANG received the Ph.D. degree in mechanical engineering from Northwestern Polytechnical University, China, in 2015. He is currently an Associate Professor with the School of Mechanical Engineering, Xi'an Polytechnic University. His main research interests include tolerance analysis, tolerance allocation, digital design, and computer-aided design. He has authored several research articles in reputed international journals. He is a Reviewer of a number of well-known scientific journals.
GUANGSHEN XU received the Ph.D. degree in mechanical engineering from Xi'an Jiaotong University, China. He is currently a Professor with the School of Mechanical Engineering, Xi'an Polytechnic University. His main research interests include tolerance analysis, tolerance allocation, digital design, computer-aided design, and 3D printing. He has authored several research articles in reputed international journals. He is a Reviewer of a number of well-known scientific journals.
SHOUJING ZHANG received the Ph.D. degree in mechanical engineering from Northwestern Polytechnical University, China, in 2009. He is currently an Associate Dean of the School of Mechanical Engineering, Xi'an Polytechnic University, China. His main research interests include manufacturing information technology and applications, intelligent manufacturing systems, industrial IoT technology and applications, smart logistics and flexible production scheduling, industrial big data technology and applications, equipment maintenance, and health management. He has been involved as a Leader and a Researcher in several national research projects in collaboration with industry. Meanwhile, he is the Senior Member of China Textile Engineering Society, and a member of China Operations Research Society.