A Data-Driven Robust Scheduling Method Integrating Particle Swarm Optimization Algorithm with Kernel-Based Estimation

The assembly job shop scheduling problem (AJSSP) widely exists in the production process of many complex products. Robust scheduling methods aim to optimize the given criteria for improving the robustness of the schedule by organizing the assembly processes under uncertainty. In this work, the uncertainty of process setup time and processing time is considered, and a framework for the robust scheduling of AJSSP using data-driven methodologies is proposed. The framework consists of obtaining the distribution information of uncertain parameters based on historical data and using a particle swarm optimization (PSO) algorithm to optimize the production schedule. Firstly, the kernel density estimation method is used to estimate the probability density function of uncertain parameters. To control the robustness of the schedule, the concept of confidence level is introduced when determining the range of uncertain parameters. Secondly, an interval scheduling method constructed using interval theory and a customized discrete PSO algorithm are used to optimize the AJSSP with assembly constraints. Several computational experiments are introduced to illustrate the proposed method, and these were proven effective in improving the performance and robustness of the schedule.


Introduction
Scheduling is one of the core issues in production management, which is very important for improving the production efficiency and resource utilization of the manufacturing system [1]. In real-world manufacturing systems, different types of uncertainties naturally appear in the production process, causing fluctuations in the production performance [2,3]. In this context, the robustness of a schedule is of great importance in addition to the schedule's quality. With the rapid development of the industrial informatics technologies, a large amount of historical data are now routinely captured and collected in many enterprises. Historical data can be used to analyze the characteristics of uncertainty, and this provides a new view for robust scheduling under uncertainty [4].
The motivation of this study is derived from scheduling problems encountered by many assembly manufacturing enterprises [5,6], where a variety of different products are obtained through multiple components by multistage machining and assembly processes (shown in Figure 1). Due to the tree-like assembly relationships, higher level components cannot be assembled until all preceding lower levels components are completed. Therefore, compared with a single machining process, fluctuations in the processing time and setup time of any part may cause the waiting of subsequent and parallel components, and this can ultimately lead to the delay of the final product [7]. In order to ensure sufficient lead time for preparing production auxiliary materials and tools, the manager should decide on the complete schedule in advance. In most real-life manufacturing systems, managers employ dispatching rules and optimization methods to generate a baseline schedule, and they then insert extra time redundancy into the schedule to ensure robustness [8,9]. However, the assessment of uncertainty in this schedule is subjective and largely depends on the experience of the manager [10]. In real-life manufacturing systems, schedules are often confronted with uncertain factors. For example, resource shortages and machine breakdowns can delay a schedule's completion time. In this work, we mainly consider the uncertainty of process setup time and processing time, which are also the two most frequent and influential uncertainties in the assembly production system. To hedge against uncertainties in the production process, traditional optimization methods, such as stochastic programming (SP), fuzzy programming (FP), and robust optimization (RO), have proposed different mathematical expressions of uncertainty [11][12][13]. The sources and types of uncertain parameters considered by these methods are different. For example, the FP method assumes that the uncertain parameters belong to a specific perturbation interval under different probabilities, and it establishes a fuzzy set based on the membership function to describe the uncertain scenarios. The SP method establishes a representative scenario tree by enumerating the possible values of uncertain parameters. These methods usually combine a priori reasoning with assumption on the uncertainty to motivate the construction of the optimization models. However, in the real world, the probability distribution of uncertainty is closely related to the specific production environment, and some complex uncertainties are even difficult to describe explicitly. Thus, traditional optimization methods are often criticized for their inflexibility or ineffectiveness in practical applications [14].
In addition to traditional optimization methods, many metaheuristic algorithms have been extensively applied to shop scheduling problems under uncertainty [15]. Among them, the particle swarm optimization algorithm is one of the most commonly used metaheuristic methods because of its flexibility and global optimization capabilities [16,17]. Han et al. [18] studied integrated production planning and scheduling with fuzzy startup time and processing time. They developed a fuzzy bi-level decision-making technique based on the PSO algorithm and a heuristic method. Jamrus et al. [19] developed a hybrid approach integrating the particle swarm optimization algorithm with genetic operators for solving a fuzzy job shop scheduling problem with an uncertain processing time. Panadero et al. [20] considered a permutation flow shop problem with stochastic processing times. They proposed a biased-randomized heuristic algorithm and extended it to a metaheuristic method by adding a variable neighborhood descent framework. Obviously, in the above research, different uncertain parameters were described and modeled by fuzzy and stochastic theories. In other studies, uncertain parameters are usually regarded as bounded uncertainties, which can be defined as a continuous interval or a finite set of discrete values. Feng et al. [21] studied a similar problem but proposed a min-max regret scheduling model with interval scenarios of processing times. They designed both exact and heuristic algorithms to solve this problem. Wang et al. [22] considered a job-shop scheduling problem with an uncertain processing time, described by discrete scenarios. Lu et al. [23] addressed a single-machine scheduling problem with uncertain processing times and sequence-dependent setup times.
All the previous studies focused on the classic flow-shop or job-shop scheduling problems under uncertainty and did not consider the assembly process. In fact, the production process with both machining and assembly is very commonly witnessed in real life. Some examples include production flows of molds, ball valves, and various complex electromechanical products [24,25]. Normally, production scheduling that considers the assembly process is more complicated, and the schedule is more sensitive to the effects of uncertainty. In addition, most of the previous studies define uncertainty based on prior assumptions or artificial experience, rather than from an objective data-driven perspective. With the recently emerging industrial informatics technologies such as the Internet of Things (IoT), sensors, and artificial intelligence, the historical data of the production process can be continuously collected and can be a new source for scheduling optimization [26,27]. These data have motivated a shift from a priori reasoning and assumptions to a new datadriven paradigm in dealing with uncertainties in the production process, which is also the focus of this research.
In the actual production process, fluctuations in the setup time and processing time of the process may cause the completion time to deviate from the expected time, and this can ultimately lead to delays in subsequent processes. Motivated by realistic needs, this paper considers two types of uncertain parameters that have an important impact on the production scheduling of the assembly manufacturing system, namely process processing time and setup time, and from there we aim to develop a data-driven robust scheduling framework. Briefly, the proposed framework consists of statistical analysis and interval scheduling methods used to describe the distribution of uncertainty parameters. On this basis, a PSO algorithm considering assembly constraints is designed to hedge against the risk of system performance degradation in uncertain environments.
The rest of this paper is organized as follows. In Section 2, we give the problem description and the notations. In Section 3, considering that the prior distribution of process setup time and processing time is unknown, we design a nonparametric estimation method based on kernel density estimation to obtain the probability density function of these two parameters. The proposed particle swarm optimization algorithm based on interval scheduling is then elaborated in detail, which is followed by a series of computational experiments and discussion. A summary covering the conclusions is presented in the last section.

Problem Description
Suppose there are several types of products with tree-like assembly structures to be produced. Each product is composed from a set of components. Those basic components produced with only machining operations are termed parts, while those assembled from other parts are termed subassemblies. The term job refers to a part or a subassembly. A final product is assembled from a set of subassemblies and parts with a sequential constraint relationship in a tree-like operation diagram. The operation tree example of a product is shown in Figure 1. The nodes in the tree operation represent the operations, while the edges represent the partial order relationship of the operation constraints. An assembly operation can only start after all required parts are completed. Each individual operation can be processed by multiple alternative machines/assembly group. For convenience, we call the above problem an assembly job shop scheduling problem.
Other assumptions made are listed as follows: • Each operation can be machined/assembled by at most one machine/assembly group at a time; • Each machine/assembly group can handle at most one operation at a time; • Only one alternative assignment is chosen for any operation; • Any started operation cannot be interrupted before it is finished; • All parts and machines are available at the beginning of the scheduling period.
In the production process, various uncertain factors need to be considered in the scheduling problem. The robustness of a schedule means that it is insensitive to various uncertainties, given a certain control policy. In this study, we assumed that managers adopt the most commonly used right-shift policy [28] to hedge against the effect of uncertainty on the baseline schedule. We used s to represent a feasible schedule that specifies the operations assigned to each machine and the corresponding sequence. Let C max (s, ξ) denote the maximum completion time of schedule s under an uncertain parameter ξ U. We considered two uncertain parameters, namely, the uncertain processing time and the setup time, and we assumed that their probability density functions can be obtained from historical data. By combining these data with the a priori structural features of uncertainties, we constructed a mathematical model for AJSSP. Some notations used for formulating the mathematical model are listed in Table 1. Table 1. Parameters and variables for modeling.

Parameters
Variables The model is presented as follows: The binary decision variable X ij is 1 if operation i is scheduled before j, and 0 otherwise. C i , which is the completion time of operation i, is ST i + SU i + PT i , the right-hand side of Equation (2). Since we are interested in minimizing the makespan, i.e., the completion time of the final operation, the objective function is max i∈I C i , which is equivalent to Equations (1) and (2). Equation (3) guarantee that each operation can only start after the related processing operations and subassembly operations are finished. Equations (4) and (5) ensure that every single machine can only process one operation simultaneously.
Obviously, the above mathematical model is a kind of semi-infinite problem and cannot be solved directly, since it has an infinite number of constraints due to the uncertain set U to which ξ belongs. Once the interval times of the schedule are determined, this model can be solved by linear programming methods, such as robust linear programming. In this research, a novel data-driven method, which combines metaheuristic algorithm and statis-tical analysis technology, is used to solve the above AJSSP under uncertainty. The details of this method are introduced in the next section.

Data-Driven Robust Scheduling Method
The proposed data-driven robust scheduling method includes two parts: an uncertain parameter analysis based on kernel density estimation technology and an interval scheduling optimization based on a particle swarm optimization algorithm. First, we estimate the probability density function of uncertain parameters based on the kernel density estimation method, and we adjust the range of the interval number through the confidence level. Second, we design an interval scheduling method that takes the interval number of uncertain parameters as input and optimizes the interval makespan through the PSO algorithm.

Kernel-Based Estimation of Uncertain Parameters
In the actual production process, uncertain parameters are usually regarded as bounded uncertainties, and managers often define the bounds of uncertain parameters based on experience. In the absence of domain-specific knowledge, the analysis of historical data provides a practical way to characterize uncertain intervals. For example, the upper and lower bounds of the interval can be specified as the maximum and minimum of data samples. However, this method may be too conservative because the worst-case scenarios at the boundary are very unlikely to occur in practice.
In this paper, we decompose the uncertain parameters SU i and PT i in the following forms: where SU i and PT represent the nominal value of the uncertain parameter;Ŝ U i and PT i denote the magnitude of the random perturbation, which can be approximated by the bounds of the samples; ξ i is a normalized random variable, which controls the range of uncertainties. By specifying the interval of ξ i , inherent uncertain parameters can be transformed to bounded uncertainties that are easy to model mathematically.
In recent years, some researchers have proposed to adopt statistical analysis and machine learning methods to learn the distributional information from the data and then establish a data-driven uncertainty set to reduce conservatism. As a well-known nonparametric estimation approach, kernel density estimation (KDE) has proven to be an effective tool for modeling probability density distribution [29]. Assume that we have the data of uncertain parameters drawn i.i.d. according to distribution P. In this paper, we propose to use the confidence region of the probability density function to quantify the uncertain interval by KDE.
Given historical data of uncertain parameters, KDE learns the probability density function of the samples according to where p(x) denotes the estimated density function of uncertain parameter x; h is a positive parameter representing the bandwidth; n is the number of samples; ϕ(·) is a kernel function, which satisfies ϕ(u) ≥ 0 (11) In our research, the Gaussian function commonly used in engineering problems is introduced as the kernel function, and the value of bandwidth h is determined through pilot experiments.
Based on the probability density functions of uncertainties, the confidence interval of the uncertain values can be calculated through interval estimation. To reduce the conservatism of the uncertain interval, we introduce an adjustable parameter σ to control the number of uncertain values that lie within the uncertain interval. Taking the uncertain parameter SU as an example, the corresponding σ can be described as represents the value ofŜ U i under the confidence level 1 − η, which can be calculated by the quartile function of SU i [30], andŜ U max i denotes the maximum range of uncertainty. When the confidence level 1 − η equals 1, all of the sampled points are contained in the uncertain interval andŜ U . Conversely, as the value of 1 − η falls to 0, fewer samples are considered (see Figure 2). Thus, the region of SU i determined by σ i can be described as: In the following optimization algorithm, the above data-driven uncertain interval is used to assist the encoding and decoding procedure.

Data-Driven PSO Algorithm with Uncertain Parameters
Over the past two decades, a large number of metaheuristic algorithms have been developed to solve production scheduling problems with complicated constraints. In this paper, we address the AJSSP with uncertain parameters, where the objectives cannot be evaluated in the form of traditional deterministic fitness functions. To solve this problem, we propose using a data-driven interval schedule to approximate objective functions under the effect of uncertainties. On this basis, this paper proposes a data-driven particle swarm optimization (DPSO) algorithm framework (see Figure 3). In DPSO, we improved the original PSO algorithm by adding the Metropolis operator and genetic operator.

Encoding Scheme of AJSSP
In the design of heuristic algorithms, the term chromosome represents the encoding of a solution. A reasonable chromosome representation is the primary task for the successful application of heuristic algorithms to solve actual production scheduling problems. In order to make the PSO algorithm applicable to the discrete AJSSP, we propose a discrete chromosome coding and initialization method based on the adjacency matrix to represent sequential constraints of the product.
To make the description about the adjacency matrix easier to follow, we assume that one product (i.e., product A) will be scheduled, and its assembly structure is shown in Figure 1. First, we establish the corresponding adjacency matrix A = a ij n×n of product A to represent the sequential constraints (see Figure 4). The row and column numbers of the matrix A represent the numbers of all parts/subassemblies to be scheduled, and A 0 is a virtual node that represents the end of assembly. In the adjacency matrix A, if the value of element a ij is not equal to 0, it means that A i is a lower level part of A j , i.e., A j can only start after all operations of A i are completed. The value of a ij represents the number of operations of part A i that have not yet been processed. In the initialization and repair process of DPSO, one can continuously search and update the adjacency matrix A to avoid generating an infeasible schedule.  Table 2 shows a feasible chromosome of product A. In the decoding process, the O-string is converted to a sequence of operations, and then each operation is assigned to the prespecified machine.

Decoding Scheme Based on Interval Schedule
By using the KDE-based estimation for historical data, the processing time and setup time of each operation can be represented as interval number PT i = PT i , PT i and SU i = SU i , SU i . The left and right bounds of the interval number can be calculated by Equation (13) and adjusted by setting different confidence levels. After that, we can construct an interval schedule based on interval theory to obtain the interval makespan.
Unlike the scheme of generating a deterministic schedule, interval theory redefines basic operation rules of interval numbers. In this study, we focus on the addition and comparison rules involved in the production scheduling process. Assuming that A and B are two interval numbers, A L /B L and A R /B R represent the left and right limits, respectively. The addition and comparison rules of A and B are defined as follows [31]: In the decoding process, Equation (15) issed to calculate the completion time of an operation, and Equation (16) is used to compare the release time of different machines.
Based on an O-string, a corresponding interval schedule can be obtained through the above operation rules. By comparing the interval numbers for C i of different operations, the interval makespan C max can be finally obtained.

Particle Update Method
In this study, we design a discrete PSO algorithm to realize the discretization of the position update process through operations such as insertion and exchange. In addition, the Metropolis criterion in the simulated annealing algorithm is used to set the acceptance probability of the location update.
To ensure the diversity of the initial population, we use a random generation method to generate initial particles. In the PSO algorithm, the particles are updated according to the interaction of particle velocity, the best-found location, and the global best-found position. In each iteration step, the position and the velocity of the i-th particle are updated in the following formulas: where pB i is the best-found location by particle i up to this step; gB i is the global best-found position among all particles up to this step; w is the inertia weight, which is used to control whether to execute h(·); RAND represents a random number between [0,1], and c is a preset positive constant. The inertia weight w is expressed in the form of probability, and w ⊗ h(X i ) means that the algorithm executes h(X i ) with the probability of w. The design of the functions h(X) and g(X, Y) refers to the commonly used genetic operations. Specifically, the function h(X) exchanges the values of two random positions in the particle X, and the function g(X, Y) is defined as the crossover operation in genetic algorithm. In addition, we introduce the Metropolis criterion to set the particle update probability in the following formulas: where X ' i represents the updated new position of particle i, and f (X i ) is the fitness function value corresponding to particle i at position X i . If i < 0, update the new position of particle i to X ' i . Otherwise, update the position of i with the probability exp(− i / f ( pB i )). In order to ensure the feasibility of the generated solutions, DPSO uses the adjacency matrix to identify whether the new solution is feasible after each operator is executed. If the solution is infeasible, the algorithm re-executes the operator until a feasible solution is obtained.

Experimental Study
To test the performance of the proposed data-driven robust scheduling method, a series of calculation experiments was performed on 2 AJSSP benchmarks and 8 cases originating from an aerospace structure parts workshop in Shanghai. Among them, Cases 1 and 2 were derived from two AJSSP examples used in Shi et al. [6]. We modified the above two cases by introducing random processing time and setup time. The features of these cases are shown in Table 3.
All the experiments were performed in Visual Studio 2017 on a personal computer with Intel Core I7-6700 2.6 GHz CPU and 4 GB RAM.

Setting Algorithm Parameters
As it is generally known, the performance of metaheuristic algorithms is greatly affected by the values of the algorithm parameters. To investigate the statistically significance of the parameter values of the algorithm, we performed an analysis of variance (ANOVA) on the obtained experimental results. In this study, 4 key parameters were considered, and each parameter contained 3 alternative value levels. Specifically, the population size included 3 levels, namely 100, 150, and 200; the inertia weight w included 3 levels, namely 0.7, 0.8, and 0.9; the constant c includes three levels, namely 0.3, 0.5, and 0.7; the function g in Equation (17) can choose three crossover operation modes, namely SPX, PMX, and OX [32]. The result of the ANOVA is shown in Table 4. In Table 4, SS means the sum of squares for the factor, which reflects the degree of variation between sample means. The abbreviation df means the degree of freedom of the factor. Here, df = k − 1, and k is the number of value levels of the factor. MS means the mean square between groups for the factor, and MS = SS/df. The label "F" is the test statistic in the ANOVA, which is equal to the ratio of the mean square between groups and the mean square within groups. The greater the value of F, the greater the difference between the groups, or in other words, the more significant the influence of this factor on the results.
According to Table 4, since the p-value (sig.) is less than the significance level (0.05), the above four parameters all have significant impact on the experimental results. Taking the relative percent difference (RPD) of the makespan as the metric, Figure 5 shows the corresponding means plot and the least significant difference (LSD) intervals of inertia weight w at the 95% confidence level for the RPD metric. When the inertia weight w is set to 0.8, the value of the RPD metric is the smallest, and thus its value is set to 0.8. Similarly, the value of population size is set to 200, while the value of constant c is set to 0.5, and the function g is set to be the SPX mode.

Robustness of Interval Schedules under Different Confidence Levels
In the proposed method, the confidence level of uncertain parameters is an important tool to control the robustness of the schedule. Based on Equation (15), different confidence levels correspond to different bounds of intervals, which make the schedule flexible to fluctuations of production parameters. Taking Case 1 as an example, the product structure and the original deterministic processing time PT i are detailed in Shi et al. [6]. In this study, the uncertain processing times PT i were assumed to follow Gaussian distribution with the mean of PT i and the standard deviation σ = 0.2 * p 0 . Here, p 0 represents the standard processing time of the operation. We obtained the optimized solution (see Figure 6) of Case 1 in the deterministic scenario and in the uncertain scenario with a 80% confidence level. As can be seen from Figure 6b, each operation in the interval schedule was represented by two intervals. Taking operation 4 as an example, the lower line O4_S indicated the interval start time of operation 4, and the upper line O4_E indicated the interval completion time, which was obtained by adding the interval processing time to the interval start time. Compared with the schedule for deterministic scenario, although the makespan of interval schedule increased, it is obvious that each operation reserves appropriate time redundancy based on its historical data, which is more suitable for the actual production process. From another perspective, the result also proves that the uncertainty of the processing time may cause the maximum 16.5% fluctuation in the makespan of Case 1, which can provide a reference for the delivery date of the product. In addition, the proposed method can help managers flexibly adjust the conservativeness and robustness of the schedules by setting different confidence levels. We tested the performance of interval schedules under different confidence levels and reported the results of Cases 1-5 in Figure 7. In Figure 7, when the confidence level was equal to 0, all parameters became deterministic values, and the interval scheduling at this time was equivalent to deterministic scheduling. With the continuous improvement of the confidence level, the interval of uncertain parameters covers more uncertain scenarios, and the value of the objective also increases. Taking Case 1 as an example, the objective value of confidence level 0 was 178, and the value of confidence level 1 was 211. This meant that the uncertainty of the processing time and setup time may cause a performance deviation of 18.5% in the makespan. In addition, for Cases 3-10, the average fluctuation of the objective value caused by different confidence levels of the processing time and setup time is 15-20%, which is basically consistent with the actual situation of the workshop.

Performance Comparison of Different Optimization Algorithms
To verify the performance of the proposed algorithm, the DPSO algorithm was compared with two classic heuristic algorithms, namely the simulated annealing algorithm (SA) and the genetic algorithm (GA) [7]. GA is a typical population-based heuristic algorithm, while SA belongs to another type of heuristic algorithm based on a single individual. We chose these two algorithms for comparison because they are the prototypes of many evolutionary algorithms of the same type, and they have shown excellent performance in many workshop scheduling problems. In addition, we also added the original PSO algorithm for comparison. The GA, SA, and PSO algorithms used the same interval-based decoding scheme as that of DPSO, and the algorithm parameters were also optimized by the ANOVA. To make a fair comparison, all algorithms used the interval numbers of uncertainty parameters under the 80% confidence interval as input, and the cut-off time of the algorithms was set to 120 s. The above three algorithms were applied on each instance with 10 independent runs. In order to comprehensively compare the stability and superiority of different algorithms, we recorded and reported the average value (denoted as Avg) and the best value (denoted as Best) of the makespan on each instance. In addition, the Avg/Best makespan of each algorithm was compared with the Avg/Best makespan of all comparison algorithms (denoted as Avg * /Best * ) from the relative percentage error (RPE). The RPE can be denoted as follows: The computational results in Table 5 show that the proposed DPSO algorithm is significantly better than the SA and GA. The pair of numbers under each case denote the scale of the case. For example, (92 × 5) means that this case contains 92 processes that need to be processed on 5 machines. For Cases 1 and 2, due to the small scale of these problems, the above three algorithms can stably obtain the same optimal solution. However, as the scale of the problem continues to increase, GA and DPSO are superior to SA in terms of optimization capability and algorithm stability. Specifically, for Cases 3 and 5, GA and DPSO can obtain the same optimal solution and perform better than SA on the same CPU time. However, even for these two cases, the average values of DPSO are 3.9% and 1.3% smaller than those of GA, which show that the stability of the proposed algorithm is better. For more complex cases (Cases 6-10), compared with the SA and GA, DPSO obtains all the optimal solutions among the three algorithms. In terms of optimization capability, the best values of optimal solutions obtained by DPSO are 3.8% and 6.3% better than GA and SA, respectively. In terms of algorithm stability, the objectives of DPSO are on average 5.1% and 10.3% better than GA and SA, respectively. Considering the similar structure of the population-based optimization algorithms, the computational results of GA and DPSO show that the proposed encoding and initialization method based on the adjacency matrix are effective. The discretized particle update method is also more advantageous than conventional genetic operations.
In addition, we compared the performance of the original PSO algorithm and the improved DPSO algorithm. For Cases 1-5, although the PSO also obtained the same best values of optimal solutions, the algorithm is not stable enough in 10 independent runs. For Cases 6-10, DPSO is obviously better than PSO in all respects. Specifically, the best values and average value of optimal solutions are increased by 2.6% and 3.3%, respectively. This result also proves the effectiveness of the proposed particle update method. In summary, the proposed DPSO algorithm is competitive in performance and stability, especially on larger-scale problems.

Conclusions
Assembly job shop is a common production organization mode in discrete manufacturing, and its production scheduling is more complicated than a single processing workshop. This paper proposes a data-driven robust scheduling method for the AJSSP with uncertain production parameters, which can generate adjustable robust schedule. The improvements of the proposed method consist of three points: (1) using the kernel density estimation method to obtain the probability density function of uncertain parameters, and introducing the confidence level to adjust the bounds of uncertain parameters; (2) establishing a decoding scheme based on interval scheduling to realize the combination of interval numbers of uncertain parameters and schedule; (3) designing a discrete PSO algorithm to deal with the assembly constraints and optimize solutions. The experimental results show that the proposed method can effectively improve the performance and robustness of the schedule. The improvement of the digital level of manufacturing enterprises is the key to applying new data analysis methods to solve traditional scheduling problems. Taking into account the complexity of the actual production environment, this study only considers the uncertainty of the processing time and setup time in the production process. Future research work can be divided into two aspects: (1) considering more types of uncertain parameters according to the actual production environment; and (2) incorporating more prior information and characteristics of uncertain parameters into the scheduling method to further improve the effectiveness and robustness of the schedule.

Conflicts of Interest:
The authors declare no conflict of interest.