Robust Fuzzy-Stochastic Programming Model and Meta-Heuristic Algorithms for Dual-Resource Constrained Flexible Job-Shop Scheduling Problem Under Machine Breakdown

Resource scheduling, job sequencing, and assigning them to available resources are the most critical issues in manufacturing systems, such as flexible job-shop systems. In addition, scheduling uncertainties have attracted significant attention in this field. This study investigates the dual-resource constrained flexible job-shop scheduling (DRCFJSS) problem under machine breakdown and operational uncertainty. Stochastic scenario-based methods were utilized to study the uncertain nature of the problem. Because process times have inherent uncertainty, they are considered fuzzy numbers and are controlled by a credibility-based measure. Robust scheduling must be developed to address unexpected disruptions, such as machine breakdowns and operational risks, such as uncertain process times. Accordingly, a novel robust fuzzy stochastic programming (RFSP) model is presented for this problem. In the proposed RFSP model, the objective function is formulated using a hybrid measure (i.e., a combined average-case and worst-case performance of the manufacturing system) under probable machine breakdown scenarios. Because the DRCFJSS problem is NP-hard, two types of meta-heuristic algorithms, evolutionary population-based, genetic algorithm (GA), and vibration damping optimization (VDO) algorithm, are used for large-sized problems. Then, the proposed RFSP model was applied to a case study, and numerical experiments with randomly generated test problems were used. In small-sized problems, the proposed model is solved using the CPLEX solver, GA and VDO algorithms. Also, the computational study confirms the proper quality of the results of the GA and VDO algorithms in medium and large-sized problems.


I. INTRODUCTION
The job-shop scheduling (JSS) problem is widely used in real-world manufacturing systems. An optimized schedule can substantially reduce costs and improve productivity. In the classical form of this problem, the best sequence for performing jobs on each machine is determined by considering the precedence and processing time constraints. It assumes that each operation of a job can be processed using a predetermined machine. However, this assumption is no longer valid for today's multipurpose machines, where a job The associate editor coordinating the review of this manuscript and approving it for publication was Muhammad Imran Tariq .
can be processed by several machines. This type of JSS problem, called the flexible job shop scheduling problem (FJSS), is more complicated than JSS, and in addition to sequencing, machine assignment problems should be considered [1].
Because the classical JSS problem has been proven to be an NP-hard problem [2], its different extensions, such as FJSS, are also considered to be NP-hard problems. Pezzella et al. [3] proposed a GA algorithm for the FJSS problem. Integrating different strategies to generate the initial population leads to better performance. Gao et al. [4] developed a hybrid GA for the FJSS problem, where a variable neighborhood descent (VND) algorithm was integrated to strengthen the search ability. Lei and Guo [5] proposed a memetic algorithm for multi-objective problems with worker flexibility, formulated as an NLIP model, aiming to minimize the maximum completion time and the individual and total maximum workloads of machines. Shen et al. [6] extended the FJSS problem (including sequence-dependent setup times) and developed a tabu search algorithm, the key elements of which were a neighborhood structure with feasibility checks and preliminary neighbor evaluations, a tabu list, and specific diversification. Kress et al. [7] also extended the problem by incorporating heterogeneous machine operator qualifications (i.e., they assumed that the processing times were machine-and operator-dependent. Some exact and heuristic decomposition-based approaches have been proposed for this problem. In fact, the scheduling problem was broken into two parts: a vehicle routing problem with precedence constraints and an operator assignment problem using logic inequalities. Liu et al. [8] proposed a hybrid distributed evolutionary model for a large-scale flexible job-shop scheduling problem (LSFJSP). The model is composed of two layers: a division layer and a coevolution layer. Three typical evolutionary algorithms were integrated into the proposed models to test their superiority.
FJSS systems that consider machine and worker constraints simultaneously, called dual-resource constrained systems, become closer to real-world situations. To optimize the performance of the system, the DRCFJSS attempts to determine the assignment of operations to machines and workers and their sequence on the assigned machines [9]. The DRCJSS problem has been considered by numerous researchers. Jaber and Neumann [10] modeled the effects of worker fatigue and recovery on the performance of a DRC system (comprising one worker). They presented a mixedinteger linear programming (MILP) model. They found that short rest breaks after each task, short cycle times, and faster recovery rates improve the system's performance. Li et al. [11] developed an adaptive hybrid approach that combines the ACO and simulated annealing (SA) algorithms. Lei and Guo [12] proposed an effective variable neighborhood search (VNS) algorithm that represents its solution as a quadruple string of ordered operations and their resources. Yazdani et al. [13] proposed an SA and VDO algorithm, both of which utilized various neighborhood structures to search the solution space. Li et al. [14] proposed a branch-population GA to minimize the makespan and cost. The branch population was incorporated to accumulate and transfer the experiences of parent chromosomes via pheromones.
Although the aforementioned extensions have made the problem more realistic, the problem becomes closer to the real-world situation by considering the uncertainties involved in the manufacturing process. The working environment may experience unexpected events such as machine breakdowns. Such events may cause the optimized schedule to stop working, and this may be infeasible under new circumstances. Furthermore, some new information, such as the due date and processing time of the new job, might be updated during the process. Accordingly, generating updated schedules that can handle the uncertainties caused by unexpected disruptions is the main concern in this field [15]. A growing body of literature addresses the uncertainties involved in manufacturing processes. Zandieh and Adibi [16] addressed a dynamic JSS problem with random job arrivals and machine breakdowns, in which a hybrid performance measure of makespan and tardiness was considered. They introduced a VNS-based algorithm and applied an event-driven policy to deal with a dynamic nature. Lei [17] minimized the makespan for scheduling a stochastic JSS problem (subjected to a breakdown). He applied the proposed GA to some test problems and compared its performance with that of the SA and PSO algorithms. Xiong et al. [18] considered two objectives for the FJSS problem, makespan, and robustness, and provided a scheduling method that is robust against machine breakdowns. This problem was solved using a multiobjective evolutionary algorithm (MOEA). Sun et al. [19] considered robustness and stability as two measures of rescheduling. They used game theory to establish an equilibrium between the two measures. Nouiri et al. [20] provided a predictive schedule that minimizes the makespan, which is robust and stable against machine breakdown. They used PSO to solve this problem. Buddala and Mahapatra [21] considered unexpected machine breakdowns for an FJSS problem. They used a two-stage teaching-learning-based optimization (2S-TLBO) method to solve the problem. Sajadi et al. [22] presented robust and stable scheduling for a DRCFJSS problem with random machine breakdowns. Two-stage GA was used to generate the predictive schedule. Meng et al. [23] addressed DRCFJSP by minimizing the energy consumption. This is the first study to study the energy consumption of DRCFJSP with a turn-off/on strategy.
The first stage of the algorithm optimized the primary objective of minimizing the makespan, where all data were assumed to be deterministic with no expected disruptions. Then, considering random machine breakdowns, the second stage optimized the makespan and stability. Tan et al. [24] proposed a fatigue-conscious DRCFJSP in their study, aiming to simultaneously relieve fatigue and promote productivity by joint scheduling of machines and workers. A multi-objective optimization model was constructed to minimize the maximum worker fatigue and makespan. Table 1 summarizes the reviewed papers. According to this table, the FJSS, and in particular the DRCFJSS problem, have rarely been addressed. In addition, these few studies on the DRCFJSS problem do not consider system disruption or machine breakdown. However, in most previous studies, the exact process time and fixed period were considered. To overcome all the above mentioned deficiencies and to fill these gaps, this research considers the DRCFJSS problem with machine breakdown and fuzzy processing time to get as close to the real world as possible.
The major contributions of this research can be summarized as follows: Dealing with machine breakdowns and fuzzy processing time, simultaneously. Scheduling both machine and worker resources in flexible job-shop manufacturing systems. Novel RFSP model for DRCFJSS problem under machine breakdown and operational uncertainty. The objective is to minimize a hybrid measure of the averageand worst-case performance of this flexible job shop system. The rest of this paper is organized as follows: In Section II, the DRCFJSS problem is described in detail. In Section III, a novel RFSP model with a hybrid objective function is presented to handle the DRCFJSS problem and cope with breakdown and uncertainty risks. Section IV describes two meta-heuristic solution approaches based on the GA and VDO to solve this NP-hard problem in large-scale instances. In Section V, a case study is presented to evaluate the proposed RFSP model and solution approach. Finally, the last section concludes the paper and presents directions for future research.

II. DRCFJSS PROBLEM UNDER MACHINE BREAKDOWN AND UNCERTAINTY
To extend the classical JSS to a DRCFJSS problem, the following definition is considered: a set of jobs (I) should be allocated to a set of machines (K) operated by a set of workers (W). Each job i consists of some operations (O i ) to be processed sequentially, which are indicated by O ij as the j th operation of job i. It is also assumed that each operation O ij can only be processed by a subset of eligible machines ⊆ K) and a subset of eligible workers (W (i,j) ⊆ L). In addition, each worker is supposed to be eligible to work on a special subset of machines (U (i,j,l) ⊆ K). Briefly, the same operation O ij may be performed by different machines, each of which is supposed to be handled by its specific certified workers. In addition, some jobs are sequence-dependent; that is, they are allowed to start immediately after their precedent ends. Regarding makespan as the schedule efficiency criteria, the DRCFJSS tries to simultaneously determine the optimal assignment of the operation to the resource (machines and workers) and optimally sequence the assigned operations of each machine.
The stochastic part of the problem originates from random machine breakdowns and uncertain job process times. Machine breakdown may occur during schedule execution and cause disruptions and unforeseen events. However, uncertain process times cause operational uncertainties during the process. It is difficult to consider the random nature of these parameters in a mathematical model. Because the main parameters have an inherent uncertainty, they are considered triangular fuzzy numbers and are controlled by a credibilitybased measure. In the proposed model, these parameters are modeled using fuzzy logic and scenario-based stochastic programming. This RFSP approach is described in detail in the following sections.
In this study, the goal is to minimize the total completion time (makespan) by considering the average-case and worst-case performances of the job-shop manufacturing system under all of the probable machine breakdown scenarios and possible values of uncertain parameters. We consider the following assumptions and constraints for the DRCFJSS problem: • All resources (machines and workers) are available at time 0.
• All jobs can be processed at time 0. • At any time, a machine can perform a maximum of one operation.
• At any time, a worker can perform a maximum of one operation.
• Each operation may be processed by a different machine and worker, but with a different processing time.
• The processing of each operation requires the worker and the machine.
• Workers may shift between their own assigned machines.
• A worker cannot shift until his current operation is completed.
• Once an operation is started on a machine, it must be completed before another operation can be started on the same machine (i.e., the preemption of the jobs is not allowed).
• Setup time and unloading time are included in the processing time.
• The transportation time between the machines is negligible or included in the processing time • Uncertain parameters such as job processing time are presented as fuzzy numbers • Machine breakdowns are presented as finite scenarios with an estimated probability of each scenario Given the latter two assumptions, a hybrid RFSP model can guarantee robust scheduling. Moreover, because the DRCFJSS problem is NP-hard, we need to develop an efficient solver based on meta-heuristic methods to solve the problem in large-scale instances. In the next two sections, the proposed robust optimization model and meta-heuristic solution approach are explained.

III. MATHEMATICAL PROGRAMMING
As mentioned above, machine breakdowns and uncertain processing times are presented as probable scenarios and fuzzy numbers, respectively. Therefore, RFSP is an appropriate approach for solving the DRCFJSS problem. We first explain the RFSP approach in general form and then present a novel mathematical model for the DRCFJSS problem using the RFSP approach.

A. ROBUST FUZZY-STOCHASTIC PROGRAMMING APPROACH WITH HYBRID MEASURE
The RFSP approach is a modern optimization under the uncertainty method used for planning under simultaneous disruption/breakdown scenarios and operational uncertainty risk [25]. RFSP is a combination of the classic scenariobased stochastic programming (SSP) approach and robust fuzzy/possibilistic programming (RPP). SSP modeling is one of the most common solution approaches utilized for scheduling, planning, and designing problems under disruption risk [26]- [28]. In the SSP approach, the decision variables are divided into two classes: scenario-independent and scenariodependent. The independent variables have the same value in all scenarios, whereas the dependent variables can be changed for each scenario. The SSP approach tends to respond to the independent variables; thus, by changing the dependent variables, the problem remains feasible and close to optimal under any probable scenario. The objective function of the SSP is to optimize the average performance of the system under all scenarios. In contrast, the RPP approach is used to cope with the inherent uncertainty and operational risk when uncertain parameters are considered as trapezoidal or triangular fuzzy numbers [29].
In general, the SSP model can be presented as follows: where z ξ is the objective value under scenario ξ ∈ , x ξ is the vector of scenario-dependent variables, and y is the vector of scenario-independent variables. Here, (c ξ , d ξ , A ξ , K ξ , b ξ ) are the parameters, (R, q) are the certain parameters, and π ξ is the probability of each scenario. If some data from optimization problems have inherent uncertainty but can be represented as a fuzzy number, fuzzy mathematical programming approaches are utilized. In fuzzy optimization, if the data are expressed as triangular fuzzy numbers or, more generally, as trapezoidal fuzzy numbers (TrFNs) (see Fig. 1), if the degree of constraint satisfaction is measured based on Possibility (Pos), Necessity (Nec), and other fuzzy measures, then a special branch called possibilistic programming is employed as follows: (4) ) is the TrFN data, and α is the confidence level based on Me (.). One of the most commonly used fuzzy measures is credibility (Cr), which is defined as the mean of two optimistic measures, Pos and pessimistic Nec. This can be calculated as follows [30]: One of the disadvantages of classical possibilistic programming is that the value of α is determined subjectively by the decision-maker and is usually adjusted passively or ultimately interactively, which may not be optimized. To overcome this, using the concepts of optimality and feasibility robustness, Pishvaee et al. [29] proposed a robust possibilistic VOLUME 9, 2021 programming (RPP) model in which the robust optimal level of α is determined proactively and objectively.
To formulate the RFSP model in general form, in the SSP model (1), b ξ is defined as a TrFN uncertain parameter (note that b ξ can also be scenario-independent). If we utilize the RPP model based on a credibility measure (3) to cope with fuzziness and the SSP approach (1) to control all probable scenarios, the RFSP model can be formulated as follows: where is a robustness factor preventing constraint violations, and ϑ (α) is total violations at the credibility level α. It is worth mentioning that ϑ (α) is defined as the difference which replaces the fuzzy number and worst-case value of TrFN (b (4)ξ ). Therefore, ϑ (α) is calculated as follows: The RFSP model (4) has a major disadvantage. It only considers the expectation value of the objective function and neglects the system performance in the worst-case scenario. To overcome this, the objective function is extended by adding the worst-case performance under all scenarios, as follows: where ρ ≥ 0 is the worst-case importance coefficient and should be set by the decision-maker. In the next subsection, we utilize the above RFSP model to formulate the DRCFJSS problem. W (C max ) Worst-case of the makespan under all of the scenarios. In the following section, we present a DRCFJSS mathematical model that considers machine breakdown and process time uncertainty. The objective function of the model is ''Min ξ ∈ π ξ .C max ξ ,'' and the constraints of the model are described (7)- (17), as shown at the bottom of the page. Equation (7) ensures that each operation is assigned to one machine and a worker. Equation (8) shows the operational sequencing of each job and calculates the completion time of the jobs. Equations (9) and (10) illustrate that if two operations from two different jobs are assigned to one machine, the processing of one of them is sequenced after the other. Similarly, Equations (11) and (12) guarantee the aforementioned sequencing rules for workers operating with each machine. Equation (13) ensures the precedence constraint for each job. Equation (14) calculates the total process time for each machine. Equation (15) calculates the number of breakdowns for each machine in each scenario. Equation (16) calculates the makespan (completion time of all jobs) under each scenario, where ξ r k .R kξ is equal to the total machine repair time. Obviously, regardless of the machine breakdowns, it is equal to the classic makespan. Equation (17) calculates the total time spent on machine repair.
To cope with the uncertainty of the above model, the proposed RFSP model for the DRCFJSS problem is presented as follows: In addition, the objective function is defined as a hybrid measure of average-and worst-case makespan minimization. (18)- (29), as shown at the bottom of page 7.

C. MODEL LINEARIZATION
Equation (18) represents the objective function of the proposed RFSP model for the DRCFJSS problem, which is the sum of the average-case makespan (M (C max )), worst-case makespan (W (C max )), and penalty for constraint violations. It is worth mentioning that and ρ are the robustness factor and worst-case importance coefficient, respectively.
Equations (21)- (27) are robust counterparts of (8)- (14) with respect to the fuzzy credibility measure at level α. According to the above formulation of the RFSP model, multiplying Y ijkl and α causes model nonlinearity. Fortunately, we can linearize this nonlinear term as follows: Equations (30) and (31) linearize the nonlinear terms of (1 − α) .Y ijkl (objective function) and (2α − 1) .Y ijkl (constraint). This linearization maintains the proposed RFSP model as a mixed linear optimization, which can be globally solved using the CPLEX solver in small-and medium-scale instances of the DRCFJSS problem. In the next section, we develop two meta-heuristic algorithms to solve the problem in large-scale instances.

IV. META-HEURISTIC SOLUTION APPROACHES
As mentioned before, the DRCFJSS problem is an NP-hard problem. Therefore, to solve the problem in largescale instances, we propose two meta-heuristic solvers: GA and VDO. GA is an evolutionary population-based meta-heuristic search algorithm that was first introduced by Booker et al. [31] and further extended by many researchers and optimization scientists. The GA is inspired by evolution phenomena using selection, crossover, and mutation Equations. (7), (15), (16) and (17) operators. In the context of FJSS problems, GA has been used in different studies (see the column of solution methodology in Table 1). In contrast, VDO is an efficient new metaheuristic local search algorithm inspired by the simulated annealing method, introduced by Mehdizadeh et al. [32] and also used by Hajipour et al. [33]. More recently, Yazdani et al. [13] indicated that the VDO solver efficiently solved the dual resource scheduling problem. We applied this method as an alternative to the GA. The efficiency of both GA and VDO solvers is highly affected by the encoding method, definition of the initial solution, and appropriate neighborhood searching operators. In the following sections, we first address the encoding structure, initial solution generators, neighborhood searching operators, and so on. Next, we illustrate two flowcharts to clarify the proposed GA and VDO solvers for the DRCFJSS problem in detail.

A. ENCODING AND SOLUTION REPRESENTATION
The chromosome structure was determined in the first step.
In some recent works, to make the search space more efficient, permutation-based approaches have been used to define the problem structure. Fig. 2 illustrates the proposed encoding for the chromosome representation. This chromosome is composed of multiple rows, each of which are chromosome genes that represent the problem variables. The first four components on the left-hand side of each row represent the selection variables, which determine the way the machines and workers are assigned to the operations (i.e., the job, operation, machine, and worker numbers, respectively). In the given example, the first row indicates that operation 3 of job 2 is assigned to machines 3 and worker 1. The fifth column was added to the matrix to determine the sequence of the operations. Each number in this column shows the number of orders in which the finishing time of the related operation is placed. In this way, a number is given to all operations while waiting in a queue. As an example, the last component in the first row indicates that operation 3 of job 2 is completed in the fourth place. Note that the chromosome length is equal to the number of operations.

B. INITIAL SOLUTION
Being valid enough in terms of quality, the generated initial solution should guarantee an acceptable level of diversity within the entire solution space. The quality of the solution is defined as its proximity to the optimal solution. To achieve the required level of these two criteria, the initial population is supposed to be generated by a combination of the following three methods:

1) RANDOM GENERATION
In this method, the initial population is generated by selecting random machines from those capable of performing any operation. Next, owing to precedence limitations, the sequence was randomly determined for the operations. This type of generation is intended to increase the level of diversity.

2) SELECTION WITH MAXIMUM PRIORITY
Zhang et al. [34] proposed a method that selects machines from a set of appropriate machines, and the operation sequence is generated based on the largest remaining total processing time. This method is intended to generate solutions with moderate diversity and quality.

3) INI-POPGEN HEURISTIC ALGORITHM
This algorithm, developed by Al-Hinai and ElMekkawy [35], is intended to generate an active schedule population. It assigns operations to machines using a determined priority based on the processing time, workload, and breakdown time of the machine. This heuristic behaves well in problems with a makespan objective function.

C. NEIGHBORHOOD SEARCHING OPERATORS
In designing the GA operators, it should be noted that producing an infeasible solution may be costly because of the run time and quality of the obtained results. For infeasible solutions, repair and penalty strategies are commonly used to return the search process to the feasible space, but both methods are time consuming and thus should be avoided. Therefore, in this study, a POX operator along with a modified mutation operator was used. In the following subsections, the implementation of each operator is described. Pursuing these steps, the search process did not exceed the feasible space.

1) POX CROSSOVER OPERATOR
Two donor parents are used as examples to describe how the operator works. The first is considered to be the donor, and the other is considered to be the receptor. First, some random genes of the job were selected from the donor parent. The selected random genes are indicated by arrows in Fig. 3. Next, the jobs of the receptor will be selected from the same job type as those selected from the donor and removed. The removed genes are shown by the circle cross marks in Fig. 3. The selected genes from the donor replace those of the receptor, and the remaining genes in the receptor remain intact. While satisfying the sequence and selection constraints, new offspring are generated; thus, infeasible chromosomes are avoided in this manner. Note that, in this way, no repair or penalty strategy is required.

2) IMPROVED POSITION BASED MUTATION OPERATOR
A mutation operator is required to generate diverse solutions. The proposed IMBM operator first selects random chromosomes as mutation candidates based on roulette wheel selection with the following probabilities: (32) where P m denotes the selection probability for the individual, and P max and P min denote the maximum and minimum mutation probabilities, respectively. The notations F, F min , and F avg denote the chromosome fitness, worst fitness among all chromosomes, and the average fitness of the population, respectively. This probability criterion allows individuals with better fitness to have a greater chance of mutating, thus extending an efficient search space that helps to avoid falling into local optima. On the other hand, individuals with lower fitness are more likely to generate better genes through crossover with better chromosomes.
After selecting the mutation candidates, a pair of genes was randomly selected from each candidate. One of the components is supposed to be removed, which is addressed by the circle cross marks in Fig. 4. The other is supposed to be reinserted into the instance owing to sequencing limitations.

D. FEASIBILITY AND REPAIR STRATEGY
As previously mentioned, repair and penalty strategies are avoided as much as possible because they are timeconsuming. However, there is still a chance to make infeasible solutions in the initial population because of the focus on diversification. Therefore, a repair strategy is adopted to check the initial solutions to regenerate feasible solutions instead of infeasible solutions. The mechanism function is called immediately after the production of each initial population. The repair mechanism determines the values of variables for randomly generated chromosomes.
Based on the constraint equations, this mechanism subtracts the two sides of each equation. The calculated value was saved as a parameter to validate the related constraints. The inequality Ax ≤ B indicates a model constraint in the general form. Subsequently, equation = B − Ax was used to calculate the subtraction value . The subtraction value must be positive to confirm its validity, and this is verified by this mechanism. If the constraint is defined as an equation instead of an inequality, then the mechanism checks whether the subtraction value is equal to zero. If any of the tests return a false answer, the chromosome must be regenerated.

E. GA AND VDO FLOWCHARTS IN SOLVING THE DRCFJSS PROBLEM
In this section, Figs. 5 and 6 show the flowchart of the proposed GA and VDO solver for solving the DRCFJSS problem.

V. COMPUTATIONAL STUDY
In this section, a case study of the DRCFJSS problem is conducted using the proposed RFSP model. In addition, the proposed model and solution approaches were evaluated by solving some experimental problems. To implement the solution methods, CPLEX 12.9 and MATLAB 2018b. software were used on a personal computer with a CPU Core TM i5 2.5 GHz and RAM 4.0 GB.

A. CASE STUDY OF THE DRCFJSS PROBLEM
First, we implemented the proposed RFSP model on a real case study related to the production line of a CNG valve factory. We studied the machining area of the factory as a DRCFJSS problem, in which different parts are assembled in the valve body. Different machine centers were used to machine the parts. The jobs that must be machined are needle nuts, needles, needle connectors, bobbin holders, pistons, magnets, nuts, fused plugs, burst discs, and excess flow nuts (N = 10). All the abovementioned parts will be processed on drilling, face machining, threading, and recessing machines according to the operations needed to complete each job (M = 6). Breakdowns may occur for any machine in the production of mechanical products. Five workers with different eligibility levels were available to operate the machines. In this manufacturing system, Table 2 lists all the jobs and operations, shows their corresponding process time as a TFN, and identifies the workers who can perform some operations. Note that the process time is considered a triangular fuzzy number, which is a special case of a TrFN. To clarify further, consider the cell at the intersection of row J3 and column O1 of this table, in which M1(2, 2.6, 2.9) represents the process time of machine M1 for operation O1 from J3. Furthermore, in this cell, W (1,2,3,4,5) identifies the workers capable of performing this operation. In Table 3, the breakdown scenarios and repair times are presented for all available machines. For example, machine M1 has two probable breakdown scenarios: one with a probability of 0.2 and two times with a probability of 0.8. It is obvious that each of the subsets O i , S (k) , V (i,j) , W (i,j) and U (i,j,l) can be derived from Table 2. In this study, none of the jobs are precedence, i.e., H (i) = ∅. The number of probable breakdown scenarios is 216 = 2 × 2 × 2 × 3 × 3 × 3. The fuzzy clustering method (FCM) can be utilized to reduce the scenarios and, consequently, the complexity problem. In this case, before implementing the proposed RFSP model, we used the FCM method to reduce the number of scenarios to the 50 most important critical scenarios.

1) NON-ROBUST SOLUTION WITHOUT USING THE RFSP MODEL
If a classical approach is used to schedule the orders of this manufacturing system, which does not consider machine breakdowns and replaces an expected/nominal value for uncertain parameters, a nonrobust solution (NRS) is obtained. The assignment of the order resources (#6 machines and #5 workers) is presented in Table 4. The sequence in which the operations are assigned to each machine is shown in Fig. 7. In the next section, we show that the NRS schedule derived from the classic/nominal approach certainly fluctuates under  different scenarios and, in each scenario, may be very different from the optimal condition. Then, the comparison of NRS with the solution obtained by the proposed RFSP model proves the superiority of the proposed robust solution.

2) ROBUST SOLUTION USING THE RFSP MODEL
In the proposed RFSP model, and ρ are two decisive factors that should be properly set to guarantee optimality robustness. Initially, parameter is set. For this, we temporarily assume that only the average performance of the manufacturing system is considered (i.e., ρ = 0). A larger results in greater robustness and decreases the possibility of constraint violation (delay). However, excessive conservatism may increase M value. Table 5 shows the M(C max ) and delay with respect to the different factors. For ≥ 16, the model robustness was 100%, and the possibility of constraint violation was 0. However, under these conservative conditions, excessive caution in allocating jobs to machines causes M(C max ) to increase significantly (26%) compared to the non-conservatism condition of = 0. By decreasing the factor, although M(C max ) decreases, the percentage of solution robustness may decrease at the same time as = 16, and the delays may increase. Therefore, ≥ 13 is an appropriate choice for the robustness factor because the robustness-or credibility-based confidence level is above 95%. Fig. 8 illustrates the sensitivity of the makespan and delay to the robustness factor, which schematically proves that = 13 is a proper setting for the robustness factor. Fig. 9 verifies that for = 13, the credibility-based robustness of the model is 95%. After adjusting the robustness factor ( = 13) in the proposed RFSP model, other influencing factors should be determined: the worst-case importance coefficient (ρ).
As mentioned earlier, the coefficient ρ in the proposed RFSP model causes W(C max ) to be considered more accurately, even if it causes W(C max ) to increase. Table 6 shows that by increasing the coefficient from ρ = 0 to ρ = 0.6, W(C max ) decreased by more than 4%, while M(C max ) increased by approximately 5%. For ρ > 0.6, although W(C max ) decreased again, the percentage increase in M(C max ) was not negligible. Therefore, ρ = 0.6 can be considered a proper choice for this factor. Fig. 10 shows the trade-off between M(C max ) and W(C max ) as a function of ρ. Fig. 11 also confirms that ρ = 0.6 is a suitable setting for the worst-case importance coefficient.
According to the above analysis, ( = 13, ρ = 0.6) is considered the best conduction of the described case study

3) EVALUATION OF THE RFSP SOLUTION
In the following paragraphs, we evaluate the robust solution obtained using the proposed RFSP model. Thus, we define two criteria to validate the RFSP model and compare it with the NSR obtained using the nominal approach. First, to explain the validation criteria, suppose that the DRCFJSS problem is solved by the ''M'' approach, and z M ξ,p is obtained from the machine breakdown under scenario ξ ∈ and the specific realization ofp for the fuzzy parameters ( p ijkl ). In addition, we define z * ξ,p as the optimum values for this breakdown scenario and fuzzy realization. Obviously, approach M is more efficient when z M ξ,p has less deviation compared to z * ξ,p and when it is feasible. In other words, lower values of z M ξ,p − z * ξ,p and p (4)ijkl −p are preferred. Therefore, we can use the following two validation criteria: where = 13 is the constraint violation penalty (robustness factor) and | | = 50 is the number of machine breakdown scenarios. The lower values of criteria 1 and 2 indicate that approach M is more efficient. The above two criteria are calculated in Table 8 for the NRS obtained by the nominal approach and the robust solution obtained by the proposed RFSP model. According to the defined criteria, the solution of the proposed RFSP approach is significantly more effective than the nominal solution TABLE 6. Sensitivity of the mean and worst makespan with respect to the worst-case importance coefficient. because it has less optimality deviation in both the average and worst cases. In fact, it improves the nominal approach solution in two criteria by 35.76% and 39.38%, respectively. Now, suppose that the parameter p ijkl is a crisp number, and we evaluate the value of the proposed stochastic approach compared to the nominal approach.  In this case, the value of the stochastic solution (VSS) is the most common and important criterion for validation. The VSS criterion helps us to specify the quality of the stochastic solution compared to the nominal/expected value. VSS is usually defined as the difference between the result of the stochastic solution and the nominal solution; the higher the amount of VSS, the more efficient the stochastic approach. Two different formulations for calculating the VSS in both the average-case and worst-case performances are as follows: where M C max * and W C max * are the mean and maximum makespan in the entire probable machine breakdown scenario, respectively, which are calculated for both the stochastic and nominal solutions. Table 9 shows these criteria and the VSS of the proposed model for different values of the worst-case importance coefficient (ρ). As shown in Table 9, the value of the second criterion (VSS WC ) is improved by an increase in ρ. However, both the VSS AC and VSS WC criteria are positive, which validates the value of the stochastic approach.

B. VALIDATION OF PROPOSED META-HEURISTIC ALGORITHMS
In the first step for the performance evolution of the proposed GA and VDO solution approach, the case study of the DRCFJSS problem is also solved using the abovementioned approaches, and the results are compared with the global optimum solution obtained by CPLEX.  Tables 10 and 11, respectively. The average S/N ratio plots for the GA and VDO algorithms are shown in Figs. 13 and 14, respectively. In addition, as a stop criterion, we consider the maximum iteration (MaxIter) limit for both the GA and VDO algorithms.
Each algorithm is performed several times, and the maximum number of iterations required for convergence is calculated and considered the MaxIter stop criterion of the algorithm.
For the above-described case study problem, the best solution of the GA is slightly different from the optimal solution obtained by CPLEX (Table 12). Fig. 15 shows the proposed GA convergence trend.
In addition to the above case study, several numerical DRFJSS test problems at different scales were randomly generated and solved by both meta-heuristic algorithms and CPLEX solvers to verify that the proposed meta-heuristic and CPLEX solvers are accurate in solving the proposed RFSP model. In addition, the maximum run-time limitation of 7200 s, as well as the maximum iteration limitation, was considered for the CPLEX solver to solve each of the experimental problems. To compare the GA and VDO with CPLEX, we first generated 10 small-sized DRCFJSS problem instances that can be optimally solved by the CPLEX solver in less than 7200 s. Then, all of the small-sized instances were solved by both the GA and VDO methods. In these small-sized instances, we calculate the optimality gap of the proposed meta-heuristics as follows: where F optimal is the optimal objective value obtained by an exact solver (CPLEX), and F metaheristics is the average result of a meta-heuristic algorithm (here GA and VDO) after multiple runs. Furthermore, we regenerate 15 medium-and large-sized instances of the DRCFJSS problem, for which the CPLEX   solver cannot find an optimal or feasible solution up to 7200 s. These problem instances were solved by the GA and VDO methods, and the relative percentage deviation (RPD) measure was applied to compare the results. The RPD was calculated as follows: where F Best is the best available value for the objective function and, as mentioned above, F metaheristics is the average value of the objective function obtained by multiple runs of a meta-heuristic algorithm. Considering the nondeterministic nature of meta-heuristic search algorithms, we perform 10 runs on the same problem instance to obtain reliable and credible results for the proposed GA and VDO. In the above PRD measure, F Best is defined as a minimum value objective function within 20 runs (10 runs of GA and 10 runs of VDO) for each problem instance. Table 13 shows that CPLEX can globally solve the proposed RFSP model for small-sized instances of DRCFJSS test problems. In these instances, the results indicate the acceptable performance of the proposed GA and VDO because the gap is negligible and the CPU run time is remarkably less than the CPLEX time, which means that the solution of the proposed meta-heuristic solvers is close to optimal, and it can be found faster.   In the medium-scale and large-scale problems, in which the CPLEX solver cannot solve optimally, each test problem is solved ten times by the GA and VDO. For each problem instance, the minimum, average, and maximum objective functions obtained by multiple runs of the proposed metaheuristics are listed in Table 14.    The RPD was calculated using (38) to evaluate the GA and VDO performances. Lower RPD values represent better performance of GA. Comparative data of the RPD measure show the improvement obtained from the GA. However, the results show that the fluctuations of the objective function values obtained by VDO are less than those obtained by the GA method, which confirms that the VDO method is more stable. Fig. 16 schematically compares the minimum, average, and maximum objective values obtained by both the GA and VDO solution approaches in 10 runs of each test problem. The objective value (M(C max )) of the GA was lower than that of the VDO.
Finally, there is no significant difference between the two proposed methods in terms of the solution time, and both methods converge in less than 7200 s. Note that the execution times of both algorithms for each experimental problem in large-sized problem instances (i.e., problems No. 11 to 25) are the same and equal to the maximum convergence time of both.
T-Student test was implemented to compare the results obtained using the GA and VDO in terms of average of objective function values. This statistical test was performed with the assumption of H0 that there was no difference in the values obtained using these two methods. As shown in Fig.17 the probability p-value of this test is higher than 0.05, indicating a decision based on the H0 hypothesis.

VI. CONCLUSION AND FUTURE RESEARCHES
In modern manufacturing systems, where a set of resources (i.e., machine and worker) processes orders/jobs, the development of a robust schedule (that optimally assigns resources to orders and determines the operation sequence) is important in two ways. First, such a schedule improves key measures, such as makespan and cost. Second, in the case of machine breakdowns or operational uncertainties, much less disruption occurs in the system performance, and as a result, it remains close to the optimal value. Thus, for the first time, the DRCFJSS problem was studied under machine breakdowns. and the inherent uncertainty associated with the processing time. To solve this problem, we assumed that breakdowns occur as a finite set of scenarios with certain chances, and the operational uncertainties were considered trapezoidal fuzzy numbers. To solve the DRCSJSS problem, a novel RFSP model was proposed, whose objective function was a hybrid measure of the average makespan and worst makespan under all scenarios. The RFSP model is formulated as a mixed linear optimization problem that can be solved using the CPLEX solver at small scales. However, because the DRCFJSS problem is NP-hard, two meta-heuristic algorithms, that is, GA and VDO, were proposed for medium and large-sized problem instances.
To validate the proposed methodology, a case study in a manufacturing system and 25 test problems of small, medium and large sizes were performed. The results verify the applicability and efficiency of the proposed model. The most important computational results are summarized as follows: • The robust solution of the proposed RFSP model is considerably advantageous to those of the classic approaches (which do not consider machine breakdowns or consider an average/nominal value for nondeterministic parameters) in terms of different measures, including average performance.
• The positivity of the VSS measure confirms the quality of the proposed stochastic solution when dealing with machine breakdowns.
• The proposed GA and VDO meta-heuristic approaches can provide a suitable solution (with a low optimality gap) for small-sized problems. In addition, they are suitable alternatives to solve medium and large-sized DRCFJSS problems.
• A comparison of the GA and VDO methods for medium and large-sized problem instances are as follows.
• From the perspective of the RPD measure, the GA performs slightly better than the VDO.
• From the perspective of stability in different runs for a given problem, the VDO performs slightly better than the GA.
• From the perspective of run time, there is no significant difference.
To develop this research and improve its applicability, it is suggested that an optimal approach be provided to reduce the number of breakdown scenarios. Moreover, we can consider the uncertainty of other parameters, such as the repair time of broken machines or, along with the makespan measure, consider other objective functions, such as the cost/profit of the manufacturing system. He is currently an Assistant Professor with the Department of Industrial Engineering, Islamic Azad University, Qazvin Branch. His research interests include production planning and scheduling, project management, applied operations research, and artificial intelligence techniques.
MAGHSOUD AMIRI received the B.S. degree from Shiraz University, Iran, in 1989, the M.S. degree from Islamic Azad University, South Tehran Branch, Iran, in 1991, and the Ph.D. degree in industrial engineering from the Sharif University of Technology, Tehran, Iran, and the Ph.D. degree in system management and productivity from Islamic Azad University, Tehran Science and Research Branch, Iran.
He is currently a Full Professor with the Department of Industrial Management, Allameh Tabataba'i University, Tehran. He is the author and publisher of 15 scientific books, an author and publisher of more than 90 articles in international journals, and 250 articles in scientific journals and conferences. His research interests include multi-criteria decision-making (MCDM), data envelopment analysis (DEA), design of experiments (DOE), response surface methodology (RSM), fuzzy MCDM, inventory control, supply chain management, simulation, and reliability engineering. He joined a group of highly cited researchers (top scientists) at the international level (ISI-ESI), in 2017, 2018, 2019, and 2020.
MOHAMMAD AMIN ADIBI was born in Iran, in 1984. He received the B.S. and M.S. degrees in industrial engineering from Qazvin Islamic Azad University (QIAU), Qazvin, Iran, in 2006 and 2010, respectively, and the Ph.D. degree in industrial engineering from the Amirkabir University of Technology, Tehran, Iran, in 2014. Since 2015, he has been an Assistant Professor at the Industrial Engineering Department, QIAU. His research interests include online optimization, optimization and data analysis interactions, dynamic scheduling, reinforcement learning, soft computing, and data analysis.