*Article* **A Multicriteria Simheuristic Approach for Solving a Stochastic Permutation Flow Shop Scheduling Problem**

**Eliana Maria Gonzalez-Neira 1,\*,†, Jairo R. Montoya-Torres <sup>2</sup> and Jose-Fernando Jimenez <sup>1</sup>**


**Abstract:** This paper proposes a hybridized simheuristic approach that couples a greedy randomized adaptive search procedure (GRASP), a Monte Carlo simulation, a Pareto archived evolution strategy (PAES), and an analytic hierarchy process (AHP), in order to solve a multicriteria stochastic permutation flow shop problem with stochastic processing times and stochastic sequence-dependent setup times. For the decisional criteria, the proposed approach considers four objective functions, including two quantitative and two qualitative criteria. While the expected value and the standard deviation of the earliness/tardiness of jobs are included in the quantitative criteria to address a robust solution in a just-in-time environment, this approach also includes a qualitative assessment of the product and customer importance in order to appraise a weighted priority for each job. An experimental design was carried out in several study instances of the flow shop problem to test the effects of the processing times and sequence-dependent setup times, obtained through lognormal and uniform probability distributions with three levels of coefficients of variation, settled as 0.3, 0.4, and 0.5. The results show that both probability distributions and coefficients of variation have a significant effect on the four decision criteria selected. In addition, the analytical hierarchical process makes it possible to choose the best sequence exhibited by the Pareto frontier that adjusts more adequately to the decision-makers' objectives.

**Keywords:** permutation flow shop; simheuristic; multicriteria; PAES; GRASP; AHP

#### **1. Introduction**

The flow shop problem (FSP) is a largely studied scheduling problem, as it models a wide set of industrial manufacturing environments [1], such as in the chemical, food-processing, automobile, and assembly industries. The purpose of solving the stated problem lies in determining the best processing sequence from n production jobs to process on m machines placed in series, aiming to optimize one or several KPIs such as the flowtime, makespan, earliness, or tardiness. A permutation flow shop scheduling problem (PFSP), which is an extension case of the FSP, considers only permutation schedules, i.e., the processing sequence of the jobs is the same for all m machines. The PFSP is considered an NP-complete problem for three or more machines [2] and an NP-hard problem for one or more machines when minimizing tardiness [3]. Furthermore, this problem augments the complexity when the uncertainty is considered for its resolution. From several approaches, a preliminary approach for tackling the stochastic permutation flow shop problems is elaborating a systematic procedure that includes uncertainties within the scheduling problems [4]. In fact, the scheduling problem and the associated complexity could be as important for designing a proper systematic technique to solve these problems [4].

After analyzing the literature on the FSP and PFSP, four aspects should be highlighted. At first, a set of approaches in the literature considers a single performance indicator as the objective function. The most common approaches are related to the makespan, due date, or

**Citation:** Gonzalez-Neira, E.M.; Montoya-Torres, J.R.; Jimenez, J.-F. A Multicriteria Simheuristic Approach for Solving a Stochastic Permutation Flow Shop Scheduling Problem. *Algorithms* **2021**, *14*, 210. https:// doi.org/10.3390/a14070210

Academic Editor: Frank Werner

Received: 27 June 2021 Accepted: 13 July 2021 Published: 14 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

just-in-time characteristics. The makespan or the expected makespan indicator, which is an absolute performance indicator, focuses on minimizing the lapsed time from the start to the end of the execution [5,6]. The due date indicator, which is a relative performance indicator, compares the objective function regarding the expected processing completion. The justin-time indicator, which is an accuracy performance indicator, indicates the positive or negative deviation occurring from an expected completion time. Certainly, the use of each performance indicator depends on the aim from the scheduling perspective.

Second, another set of approaches is focused on single-objective problems rather than multicriteria problems. Certainly, even industrial problems must solve various objectives simultaneously [7], and researchers have focused on optimizing a single objective function rather than multiple ones. This issue can be addressed by considering earliness and tardiness objectives jointly in a JIT environment.

Third, most studies consider only quantitative decision criteria. However, researchers have recently become interested in qualitative criteria as this approach might reduce the gap between the theoretical concept of problems and the execution. Some examples are Chang and Lo [8] and Chang et al. [9]. These authors studied a multicriteria job shop, whereas the customers' strategic importance was considered as a qualitative criterion. While Chang and Lo [8] proposed a hybridized genetic algorithm, tabu search, analytic hierarchy process, and fuzzy theory to solve the problem, Chang et al. [9] proposed a hybridization of an ant colony algorithm and an analytic hierarchy process for solving the FSP. Another approach minimized the expected costs of tardiness as a quantitative criterion and strategic customer importance as a qualitative criterion in a stochastic hybrid FSP [10]. The authors employed a GRASP metaheuristic and a Monte Carlo simulation method with a stochastic multicriteria acceptability analysis to handle both qualitative and quantitative criteria.

Finally, some research approaches consider scheduling under uncertain conditions, and these are generally divided into two approaches: the stochastic approach, in which parameters are modeled with probability distributions (PDs) aiming to minimize the expected value of a selected metric, and the robust approach (RA), in which uncertain parameters are modeled with intervals and the schedule obtained is more stable and less variable. Nonetheless, previous studies have not dealt with a combination of stochastic and robust approaches for solving the flow shop problem. Companies can collect production data in a short time, yielding enough data to accurately estimate the probability distribution of uncertain parameters. Then, we believe that the latter approach leverages the robust schedule as it can be more easily adjusted than other schedules in which uncertainties are modeled with intervals.

One of the recent approaches to solve stochastic combinatorial optimization problems is simheuristics. These hybridize a metaheuristic approach with a simulation to obtain solutions for stochastic problems. Simheuristics have been successfully used in vehicle routing problems [11–13], inventory routing problems [14], facility location problems [15], and scheduling problems [16–19].

In this sense, this paper attempts to contribute to the literature by proposing a systematic technique for solving the stochastic permutation flow shop problem (SPFSP) considering stochastic processing times and sequence-dependent setup times to optimize multiple criteria. To the best of our knowledge, no previous study has included the simultaneous analysis of a JIT environment with stochastic parameters, quantitative metrics, and qualitative metrics for obtaining robust solutions. Therefore, the current research proposes a multicriteria simheuristic approach to solve an SPFSP, including both quantitative and qualitative objectives, providing robust solutions. On the one hand, for quantitative objectives, the proposed approach addresses the expected earliness/tardiness (E[E/T]) and the standard deviation of earliness/tardiness (SD(E/T)), for estimating the JIT metrics and obtaining several robust schedules. On the other hand, as qualitative metrics, this research considers the expected customer importance of jobs (E[CI]) and the expected product importance (E[PI]), favoring the job priority for the company. Then, this paper considers both

customer and product importance as a company might be more interested in delivering high-profitability products regardless of the customer or, conversely, it might be more interested in delivering frequently with high-spending customers rather than sporadic ones. This paper is an extension of a previous conference work [20] that only considered one qualitative criterion and solved a partial set of benchmark instances. This work includes one more qualitative criterion (the maximization of the accomplishment of product importance) and executes more computational experiments by evaluating 180 instances proposed by [21] and analyzing the behavior of Pareto frontiers with more coefficients of variation of processing times and sequence-dependent setup times.

The remainder of this paper is organized as follows. Section 2 presents a literature review of the single-objective stochastic FSP (SFSP), robust FSP, and multi-objective SFSP. Section 3 describes the proposed solution for the SPFSP under qualitative and quantitative decision criteria. Section 4 provides the results of the computational experiments that validate the proposed approach. Finally, conclusions and recommendations for future work are presented in Section 5.

#### **2. Literature Review**

The deterministic FSP is one of the most frequently studied problems in the scheduling literature, but the stochastic version has been addressed less often. The FSP under the conditions of uncertainty has received more attention recently because it is closer to real manufacturing environments. Nevertheless, there are fewer FSP studies under uncertain conditions than there are deterministic FSP studies. For the literature reviews, deterministic FSPs and its solution methods have more than ten literature reviews, including Yenisey and Yagmahan [7], Pan and Ruiz [22], Arora and Agarwal [23], Nagano and Miyata [24], Rossit et al. [25], and Fernandez-Viagas et al. [26]. Meanwhile, their stochastic and uncertain counterparts have had only two literature reviews in the past 18 years, Gourgand et al. [5] and González-Neira et al. [27]. Studies of the FSP with uncertainty generally use one of three approaches to deal with uncertain parameters: the stochastic approach, the robustness approach, and the fuzzy approach. The most frequently used approach is stochastic, with the uncertain parameters modeled using probability distributions. Examples of this approach can be found in Framinan and Perez-Gonzalez [28], Lin and Chen [29], and Qin et al. [30]. For the robustness approach, there is no need for acknowledging the distribution of the data of uncertain parameters due to this being modeled with intervals or datasets. Examples of this approach can be found in Fazayeli et al. [31] and Ying [32]. For the fuzzy approach, uncertain parameters are modeled using fuzzy numbers, as in Behnamian and Fatemi Ghomi [33] and Huang et al. [34].

The aim of this paper is to obtain robust schedules for the stochastic permutation flow shop problem (SPFSP), considering both quantitative and qualitative objective criteria. Then, the literature review in this paper is organized into four subsections: Section 2.1 presents the literature review for single-objective SFSPs. Section 2.2 presents an overview of studies using the robust FSP with a single objective. Section 2.3 examines current approaches for a multi-objective FSP under uncertain conditions, including robust, stochastic, and fuzzy approaches. Section 2.4 presents a review of qualitative production criteria.

#### *2.1. Single-Objective SFSP*

Since the 1950s, the deterministic FSP has received considerable attention, but the SFSP has been less studied because it is more difficult to solve than the deterministic version. Nevertheless, with the advances in computation, the SFSP has been studied more in the last few years. Therefore, we present in chronological order the most important studies on the SFSP since its beginnings.

From the early days of this topic until the 2010s, the following approaches have been found. Makino [35] presented an optimal solution for the expected makespan in an SFSP with two jobs and general distributions for two machines, as well as with three machines, considering the exponential and Erlang distributions of the processing times. Talwar [36] proposed a rule to optimally sequence n jobs and two machines with exponentially distributed processing times, and this was proved optimal by Cunningham and Dutta [37]. Alcaide et al. [38] developed a dynamic procedure for considering stochastic machine breakdowns in an SFSP. This procedure can obtain an optimal solution for the initial stochastic problem when the associated without-breakdowns stochastic partial problems are solved optimally. Gourgand et al. [39] proposed simulated annealing with a Markovian model to solve the SPFSP with stochastic processing times and limited buffers. Wang et al. [40] presented a genetic ordinal optimization approach that hybridizes ordinal optimization, optimal computing budget allocation, and a genetic algorithm with uniformly distributed processing times. Kalczynski and Kamburowski [41] developed a new rule that enables optimal solutions when processing times are Weibull distributed. This rule includes Johnson's and Talwar's rules as special cases. Portougal and Trietsch [42] proposed a heuristic called API that consists of two steps: (1) the ordering of jobs with expected processing times through Johnson's rule and (2) applying all possible interchanges of two jobs in the sequence.

Then, in the early 2010s, Baker and Trietsch [43] compared Talwar's, Johnson's, and the API rules, testing different types of probability distributions. The authors concluded that for a high coefficient of variation, Talwar's and Johnson's rules yield better results, but when the coefficients of variation are low, Johnson's rule alone provides good results. Baker and Altheimer [44] compared three heuristic procedures adapted from rules that performed well in deterministic counterparts. The procedures used were CDS/Johnson, CDS/Talwar, and NEH. The authors found that NEH had the best results. Elyasi and Salmasi [45] solved the SFSP considering normally distributed processing times and variances proportional to their means, as well as gamma distributed processing times with the same scale parameter for all jobs. The authors aimed to minimize the expected tardy jobs by using chance constraint programming. Elyasi and Salmasi [46] proposed a dynamic method to solve the SFSP in which the stochastic parameters were the due dates of jobs, in order to minimize the expected tardy jobs. Juan et al. [18] developed a simheuristic that involves an iterated local search metaheuristic considering stochastic processing times distributed lognormally. Framinan and Perez-Gonzalez [28] compared well-known heuristics through a simulation procedure with a variable number of iterations and different probability distributions. The procedures compared were stochastic NEH, stochastic CDS/Talwar, deterministic NEH, deterministic CDS/Talwar, and deterministic NEH using CDS/Talwar as the initial solution. These authors found that stochastic NEH obtained the best results.

Finally, in the late 2010s and the beginning of the 2020s, Gonzalez-Neira et al. [47] minimized the expected makespan in a distributed assembly PFSP with stochastic processing times through a GRASP simheuristic, evaluating the robustness of the problem as well. González-Neira et al. [48] compared thirteen dispatching rules through a simulation procedure with variable iterations for ten different objective functions (five expected values and five standard deviations of the makespan, flowtime, tardiness, maximum tardiness, and tardy jobs) under lognormal, uniform, and exponential distributions of processing times. Hatami et al. [17] addressed a parallel SFSP where products had components produced in different parallel flow shops. The authors developed a simheuristic that embeds an iterated local search algorithm to minimize the expected makespan and makespan percentiles. Marichelvam and Geetha [49] considered uncertain processing times and machine breakdowns to minimize the makespan. To solve the problem, a hybridization of a firefly algorithm and a variable neighborhood algorithm was designed, demonstrating promising results with extensive computational experiments. Villarinho et al. [50] proposed a simheuristic that integrates a biased randomized heuristic into a variable neighborhood descent with Monte Carlo simulation to maximize the expected payoff in an SFSP that considers stochastic processing times.

#### *2.2. Robust FSP*

Such as the studies on the SFSP, few studies have addressed the robustness for the FSP. Table 1 presents the main characteristics of previous studies in this field. As the table shows, most of them studied the makespan as a single objective, and only three works addressed multiple objectives. These three studies included only one parameter under uncertainty conditions. While Liu et al. [51] studied stochastic machine breakdowns, dynamic job arrivals, and unexpected job availability, Liao and Fu [52] and Goli et al. [53] considered uncertain processing times. It is important to note that none of these works simultaneously studied stochastic sequence-dependent setup times and stochastic processing times while considering robustness.

#### *2.3. Multi-Objective FSP under Uncertainty Conditions*

Few works have examined uncertainties simultaneously with multi-objective decisions in comparison with a single objective. Table 2 presents the main characteristics of the studies performed in this area. Makespan, tardiness, and flowtime are the most analyzed measures, and earliness is not widely studied. In addition, stochastic setup times are not considered.


**Table 1.** Modeling approaches for the stochastic flow shop problem (SFSP) that address robustness.


good results in comparison to another algorithm.

**Table 2.** Modeling approaches for the multi-objective stochastic flow shop problem (MO-SFSP).

#### *2.4. Qualitative Criteria*

Within organizations, main processes have their own objectives, which are often in conflict. For instance, the objectives of marketing involve maximizing service levels and sales; procurement seeks the prioritization of products and order replenishment; and production and manufacturing aim to maximize throughput and minimize costs. That is why it is important to use an approach that takes into account all of these objectives simultaneously [72]. In this context, production should consider marketing criteria (often qualitative) in the decision-making process in order to improve customer service and reduce conflicts between marketing and production [73]. Among the marketing objectives, production planning is affected by the product importance and customer importance. Georgakopoulos and Mihiotis [74] analyzed these two aspects in a distribution network design. Considering the product importance, aspects such as turnover, profit rate, image, and discount policies must be considered in order to categorize the product. With regard to customer importance, aspects such as turnover, image, and customer requirements must be considered in order to categorize the customer. For instance, González-Neira et al. [10] included the customer importance in a hybrid FSP through the integral analysis method García Cáceres et al. [75] based on stochastic multicriteria acceptability analysis Lahdelma et al. [76].

The academic literature includes very few studies on scheduling problems that considered qualitative decision criteria. As stated above, Chang and Lo [8] and Chang et al. [9] proposed a multicriteria objective function that included qualitative aspects such as marketing considerations, the strategic importance of customers, and order profit/risk in a job shop environment with fuzzy parameters. Chang and Lo [8] proposed a GA/TS approach hybridized with AHP, while Chang et al. [9] examined an ant colony optimization with an AHP process to solve job shop problems. In our opinion, some possible reasons for the few studies on this research topic are the lack of objectivity while measuring these metrics and the difficulty of having unbiased indicators.

#### **3. Proposed Approach**

This paper proposes a simheuristic technique that integrates Monte Carlo simulation into a GRASP metaheuristic hybridized with the Pareto archived evolution strategy (PAES) technique for optimizing multiple objectives. Interested readers may consult Resende and Ribeiro [77] for the GRASP metaheuristic and Knowles and Corne [78] for the PAES technique. In addition, an AHP methodology is integrated with the simheuristic to assess the Pareto solutions under different weight arrays for the selected criteria. The algorithm is named multicriteria simheuristic with GRASP (MC-SIM-GRASP). A GRASP metaheuristic for optimization purposes was selected due to it having the advantage of constructing its own initial solution and the no memory characteristic, which makes it useful for scheduling problems. PAES was selected because it has the advantage of avoiding favoring a search direction in the local search and exploring different solutions across the Pareto front.

#### *3.1. MC-SIM-GRASP: Construction Phase*

The purpose of the construction phase in GRASP and in the proposed algorithm is to construct a solution by the sequential aggregation of jobs to the solution from a set of possible jobs ranked through a greedy or fitness function. Even though the studied problem is a multicriteria objective, a pure strategy was selected for the construction phase, meaning that a unique objective function guides the entire construction [79]. Then, for this phase, an earliest due date (EDD) rule was selected to deal with the earliness/tardiness objective as a single greedy function. However, considering the other functions, a respective penalization was included for the customer and product importance functions regarding each job, each depending on customer importance or product importance, in comparison with the position of the job in the sequence. Table 3 shows the penalization for the customer importance criterion for an example of 10 jobs with five levels of customer importance. These penalization scores are based on the following criteria: A job that is not belated has a score of zero. If a job is delayed, the penalization is greater if the customer importance

is high and lower if the customer importance is low. Then, the job penalization increases if the job is taking the place of a job that has greater importance. The case study for this paper defines five levels of customer importance, where Level 1 corresponds to the most important customers and Level 5 to the least important. For the instances tested in this research, the customer importance for each job was randomly assigned using the probabilities indicated in Table 4. Naturally, this scale for customer importance and the probability of the importance level were established here for testing purposes. In real scenarios, the assignment of customer importance will not be probabilistic, but rather deterministic according to the views of the decision-maker. A similar procedure was conducted to assign a product importance for each job and the corresponding penalization.


**Table 3.** Job sequence position penalization depending on customer importance.

**Table 4.** Customer importance probability.


The main contribution of the construction phase of the proposed approach is the alternation of three different greedy functions at each iteration of GRASP. For instance, the first iteration begins the construction phase with an EDD rule criterion; the second iteration continues using the customer importance; and the third iteration continues using the product importance. These iterations are repeated in the same order until the solution/schedule is completed. It is important to note that the EDD rule criterion is used to guide the constructions for both the earliness/tardiness mean and standard deviation objectives at the same time. Inside each step of construction, the restricted candidate list (RCL) is defined as the subset of jobs with the best value of 10% of the total range of the greedy function values obtained. A job is then randomly selected from the RCL to form part of the solution. The method continues iteratively until all jobs have been scheduled, included in the solution. After the sequence is completed, a Monte Carlo simulation is performed with as many runs as needed to obtain confidence intervals of at least 1% for the four objective functions (expected earliness/tardiness, standard deviation of earliness/tardiness, customer importance, and product importance). In the first iteration of the construction phase, this solution is included in the nondominated solutions (NDS) archive, which is initially empty. In the remaining iterations, this solution is evaluated using PAES to determine whether it is considered within the set of NDS. Figure 1 presents the flow diagram of the construction phase of the proposed approach.

**Figure 1.** Flowchart of the GRASP construction phase.

#### *3.2. MC-SIM-GRASP: Local Search Phase*

The purpose of the local search phase in GRASP and in the proposed algorithm is to execute consecutive iterative improvements of the solution obtained from the construction phase to obtain a better solution. In this paper, the local search phase consists of 2 opt interchanges between jobs. Each time an interchange is performed, a Monte Carlo simulation is executed to estimate the four objective functions with a confidence interval accuracy of at least ±1% for each value and a confidence of 95%. The solution is then evaluated to decide whether it should be placed in the NDS archive. If so, the other solutions already saved in the NDS archive are evaluated to determine whether they should remain in the NDS archive. If the solution does not belong in the NDS archive, it is discarded. This procedure is conducted according to the PAES proposed by Knowles and Corne [78]. A MC-SIM-GRASP iteration ends when no interchanges can enter into the NDS archive, and a new iteration is begun, maintaining the actual NDS archive. The simheuristic stop time is established as the number of jobs × number of machines × 1.0 second. Figure 2 presents the flowchart of the local search phase.

**Figure 2.** Flowchart of the GRASP local search phase.

#### *3.3. NDS Archive Solution Selection Using the AHP Methodology*

After the MC-SIM-GRASP has finished, the entire Pareto sequence is scored using the AHP methodology. As is shown in Table 5, eight vectors are used for the criteria weights for the four objectives. First, a 4 × 4 pairwise comparison matrix is created for scoring each *i*th objective function versus the *j*th objective, on a scale from 1 to 9, as indicated in the AHP technique. An example of the obtainment of a criteria weight vector is shown in Table 6. Then, eight different comparisons are performed to obtain the eight mentioned different vectors of the criteria weights. Considering the usage of each weight vector, the Pareto frontier solution that presents the best AHP score can be selected. In order to compute the matrix of option scores, for each pair of sequences *s*<sup>1</sup> and *s*2, we divided the expected earliness/tardiness of *s*<sup>1</sup> by the expected earliness/tardiness of *s*2, so if the division is greater than 1, the earliness/tardiness of *s*<sup>1</sup> is worse than the earliness/tardiness of *s*2, and vice versa. Similar divisions are performed for the other two objective functions (standard deviation of earliness/tardiness and customer importance).

**Table 5.** Criteria weights for Pareto solutions.



**Table 6.** Example of AHP scores and the resultant priority vector.

#### **4. Computational Experiments and Statistical Analysis**

For experimentation purposes, a set of 180 benchmark instances proposed by Ciavotta et al. [21] were selected to evaluate the effects of different probability distributions (PDs) and coefficients of variation (CVs) of the processing and setup times in the objective functions. Specifically, the 180 instances were 10 instances for each combination of 20, 50, and 100 jobs with 5, 10, and 20 machines and a size of sequence-dependent setup times of 50% and 125%. Two PDs, lognormal, uniform, and three CVs, 0.3, 0.4, and 0.5, were selected to model both the stochastic processing and setup times. This corresponds to 6480 NDS archives. The proposed method was implemented in Java and run on an Intel Core i7-4770 with a 3.4 GHz processor and 8GB of RAM. The stopping criterion for MC-SIM-GRASP was established as *numberO f Jobs* · *numberO f Machines* · 1.0 s. The best-qualified solution of each NDS archive was selected by using the AHP method. This results in a total of 51,840 solutions, each of which has an AHP score and a value for the four objective functions. Four ANOVAs were conducted to jointly analyze the effect of the eight factors on the four selected objective functions (E[E/T], SD[E/T], E[CI], and E[PI]). The factors selected for the experimental design are: probability distribution of processing times (PDPT), coefficient of variation of processing times (CVPT), probability distribution of setup times (PDST), coefficient of variation of setup times (CVST), the vectors of criteria weights of the AHP (WV), number of jobs, number of machines, and generation size of the standard deviation of sequence-dependent setup times (SDST). The factors and their levels are presented in Table 7.


**Table 7.** Factors of the experimental design and their corresponding levels.

The results showed that all main effects are statistically significant on the four objective functions except PDPT for E[CI] and E[PI] and the PDST and CVST for E[PI]. Nevertheless, some double and triple interactions that include the PDPT, PDST, and CVST have a significant effect on E[CI] and E[PI]. In fact, for at least one of the four objective functions, the double interaction effects are also significant (see Table 8). This shows that the WV discriminates among the Pareto solutions, helping the decision-maker select a solution from the Pareto frontier. We identified in addition that as the CVPT and CVST increase, the expected value of earliness (E[E/T]) and the standard deviation of earliness tardiness indicator (SD[E/T]) augment as well. The same occurs with E[CI], but not to the same degree. Additionally, the measures tend to be greater for the lognormal distribution than for the uniform distribution. This shows the importance of accurately fitting the PD to obtain adjusted robust measures.

Figures 3 and 4 show the main effect plots of the factors WV, PDST, PDPT, CVST, and CVPT on E[E/T] and SD[E/T], respectively. The axes of the main effect plots are the levels of each factor. It can be seen that for different weight vectors, the objectives E[E/T] and SD[E/T] imply that the AHP is capable of selecting a different solution of the Pareto frontier according to the preferences given by the decision-maker in the AHP method. In the case of the probability distributions (the PDST and PDPT factors), the lognormal distribution presents a slightly greater E [E/T] and SD[E/T] in comparison with the uniform distribution, which gives the idea that despite using the same expected values of the processing and setup times under the same coefficient of variation, the solutions vary with the change of the distribution used. Additionally, as expected, as the coefficient of variation of the setup and the processing times increase (the CVST and CVPT factors) the objectives E[E/T] and SD[E/T] also increment, indicating a higher variability in the obtained solutions of the Pareto frontier. These aspects led us to highlight the importance of including uncertainty in the optimization problem when it is really present, and the value of making an accurate distribution fitting to make adequate decisions.


**Table 8.** *p*-values and *R*<sup>2</sup> *adj* of the ANOVAs for each objective function.

**Figure 3.** Main effects plots for E[E/T].

**Figure 4.** Main effects plots for S(E/T).

Figures 5 and 6 present the interaction plots between factors CVST and CVPT for the E[CI] and E[PI] objectives, respectively. The axes of the plots are the levels of the coefficients of variation. These plots confirm that for both qualitative objectives, there exists an interaction effect between the coefficients of variation of the processing times (CVPT) and the coefficients of variation of the setup times (CVST), which means that the best solution selected with the AHP method, in terms of the qualitative criteria, varies depending on the variability of the processing and setup times. This leads again to the relevance of including uncertainties in the optimization process and the execution of an accurate distribution fitting.

**Figure 5.** Interaction plots of the CVST-CVPT for E[E/T].

**Figure 6.** Interaction plots of the CVST-CVPT for S(E/T).

To conduct the experimental design and validate the validity and objectivity, the assumptions of homoscedasticity, normality, and independence were tested. Because the homoscedasticity and normality tests were not fulfilled, we performed a Friedman test to corroborate the ANOVA results. Tables 9 and 10 present the detailed Friedman tests for E[E/T] and SD[E/T] in terms of the factors PDPT, CVPT, PDST, and CVST. To the best of our knowledge, this is the only work that has studied these four objective functions in an SPFSP. We present three indicators for the multi-objective problems proposed by Ebrahimi et al. [80] and Karimi et al. [81], which can be used for future comparisons. For this work, these indicators were adjusted for the four objective functions analyzed:


$$MID = \frac{\sum\_{i=1}^{NPS} c\_i}{NPS} \tag{1}$$

where:

$$c\_i = \sqrt{E[E/T]\_i^2 + SD[E/T]\_i^2 + E[CI]\_i^2 + E[PI]\_i^2} \tag{2}$$

$$SNS = \sqrt{\frac{\sum\_{i=1}^{NPS} (MID - c\_i)^2}{NPS}} \tag{3}$$



**Table 10.** Friedman test for SD(E/T) versus the factors PDPT, PDST, CVPT, and CVST.


Tables 11–13 show the averages of the NPS, MID, and SNS for each instance size, each combination of the PDPT with CVPT, and each combination of the PDST with CVST, respectively. In Table 11, it can be seen that the MID and SNS increase as the number of jobs or machines increases, and due to the increment of the jobs or machines, the expected and standard deviation of tardiness present a crescent behavior. Moreover, the NPS also tends to augment, giving more possible solutions to chose in scenarios with higher variability.


**Table 11.** Performance results of the proposed approach.

Table 12 shows the minimum difference in the MID and SNS between the lognormal and uniform distributions for the processing times, under the same coefficient of variation. Nevertheless, as the coefficient of variation of the processing times increases, independently of the probability distribution, the MID and SNS augment. That means that the quality of the Pareto frontier is worse as the coefficient of variation of the processing times increases. This suggests that production managers should encourage a continuous improvement of the production processes to reduce the variability of the process insofar as that is possible. Moreover, it is interesting to see that the NPS decreases as the coefficient of variation of the processing times increments, showing that high variability scenarios present fewer possible solutions to choose for the decision-maker.


**Table 12.** Performance multi-objective metrics for the PDPT and CVPT.

The measures presented in Table 13 show a large difference in the SNS between the lognormal and uniform distributions of the sequence-dependent setup times, under the same coefficient of variation of the setup times, which is a different behavior than that presented for the coefficient of variation of the processing times. Moreover, the NPS also presents high differences between the uniform and lognormal distributions of the setup times, whereas for the coefficient of variation of the processing times, the NPS values were very close. This allowed us to conclude that each input parameter can cause different performances of the solutions, and thus, the AHP selects the corresponding solution of the Pareto frontier to fulfill the decision-maker's expectations.


**Table 13.** Performance multi-objective metrics for the PDST and CVST.

#### **5. Conclusions and Recommendations**

This paper presented a multicriteria simheuristic that hybridizes a GRASP, a PAES algorithm, and a Monte Carlo simulation (MC-SIM-GRASP) to solve a multi-objective stochastic permutation flow shop scheduling problem (SPFSP). The approach obtains a set of nondominated solutions for four objectives: for expected earliness/tardiness, standard deviation of earliness/tardiness, expected customer importance, and expected product importance in an SPFSP. It was used an analytical hierarchical process (AHP) to select the desired solution for the decision-maker from the set of solutions from the Pareto frontier. The purpose of this method was to ease the selection of a solution and include a qualitative criterion for the objective of the decision-making. This paper analyzed the effect of eight factors in the behavior of the four objective functions selected. To this end, four ANOVAs were carried out with seven factors: type of probability distribution (PD) of stochastic processing times, coefficient of variation (CV) of stochastic processing times, PD of sequence-dependent setup times, CV of sequence-dependent setup times, the vector weights of criteria in the AHP methodology, the number of jobs, and the number of machines. The outcomes showed that all factors had significant effects on the four objective functions. The results obtained in this paper support the importance of including uncertainty, modeled with adequate PDs, to obtain robust solutions. Additionally, a set of multi-objective metrics (number of Pareto solutions, means' ideal distance, and spread of the nondominance solution) was calculated for future comparisons, as this problem has not been solved in the literature before. Future work may analyze other PDs and CVs. It would be useful to analyze a case in which the processing time PD of each job has a different CV, which is generally true in real cases. Finally, other qualitative criteria should be incorporated in the analysis.

**Author Contributions:** Conceptualization, E.M.G.-N.; methodology, E.M.G.-N.; programming, E.M.G.-N.; formal analysis, E.M.G.-N.; writing—original draft, E.M.G.-N. and J.-F.J.; writing—review and editing, J.-F.J. and J.-F.J.; supervision, J.R.M.-T. All authors read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The benchmark instances evaluated in this paper were taken from Ciavotta et al. [21], and the results obtained for the analysis are available in https://livejaverianaedu -my.sharepoint.com/:x:/g/personal/eliana\_gonzalez\_javeriana\_edu\_co/EdOPPW3UqvtLsAmLuso ECKwBZXGoMPwwfHQFQCCkbh1-Qg?e=54af92 (acessed on 27 June 2021).

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

#### **References**

