Minimizing the total tardiness and makespan in an open shop scheduling problem with sequence-dependent setup times

Noori-Darvish, Samaneh; Tavakkoli-Moghaddam, Reza Article Minimizing the total tardiness and makespan in an open shop scheduling problem with sequence-dependent setup times Journal of Industrial Engineering International Provided in Cooperation with: Islamic Azad University (IAU), Tehran Suggested Citation: Noori-Darvish, Samaneh; Tavakkoli-Moghaddam, Reza (2012) : Minimizing the total tardiness and makespan in an open shop scheduling problem with sequence-dependent setup times, Journal of Industrial Engineering International, ISSN 2251-712X, Springer, Heidelberg, Vol. 8, pp. 1-13, https://doi.org/10.1186/2251-712X-8-25 This Version is available at: https://hdl.handle.net/10419/78593 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. http://creativecommons.org/licenses/by/2.0/ ORIGINAL RESEARCH Open Access Minimizing the total tardiness and makespan in an open shop scheduling problem with sequence-dependent setup times Samaneh Noori-Darvish 1 and Reza Tavakkoli-Moghaddam 2* Abstract We consider an open shop scheduling problem with setup and processing times separately such that not only the setup times are dependent on the machines, but also they are dependent on the sequence of jobs that should be processed on a machine. A novel bi-objective mathematical programming is designed in order to minimize the total tardiness and the makespan. Among several multi-objective decision making (MODM) methods, an interactive one, called the TH method is applied for solving small-sized instances optimally and obtaining Pareto-optimal solutions by the Lingo software. To achieve Pareto-optimal sets for medium to large-sized problems, an improved non-dominated sorting genetic algorithm II (NSGA-II) is presented that consists of a heuristic method for obtaining a good initial population. In addition, by using the design of experiments (DOE), the efficiency of the proposed improved NSGA-II is compared with the efficiency of a well-known multi-objective genetic algorithm, namely SPEAII. Finally, the performance of the improved NSGA-II is examined in a comparison with the performance of the traditional NSGA-II. Keywords: Open shop scheduling, Total tardiness, Makespan, Sequence-dependent setup times, NSGA-II, SPEA-II Background An open shop scheduling problem (OSSP) is a kind of shop scheduling such that the operations can be executed in any order. The open shop allows much flexibility in scheduling, but it is difficult to develop rules that give an optimum sequence for every problem (Sule 1997). This problem is a class of NP-hard ones (Gonzalez and Sahni 1976. In this paper, we consider a special feature in OSSPs, called sequence-dependent setup time. The process of preparing machines between jobs is considered as a setup. In fact, setup times affect on the completion time of each job. As a result, they also affect on tardiness, earliness and other important criteria. Allahverdi et al. (2008) surveyed the literature of setup times or costs in scheduling problems. They classified scheduling problems into those with batching and nonbatching considerations as well as sequence-independent and sequence-dependent setup times. According to the technology and the kind of machines used in a work environment and variety of products, setup times can be dependent on both machines and the sequence of jobs that should be processed on a machine. In many practical production systems (e.g., chemical, printing, pharmaceutical and automobile manufacturing), the setup tasks (i.e., cleaning up and changing tools) are sequence-dependent (Zandieh et al., 2006 and Roshanaei et al., 2009). Low and Yeh (2008) addressed an open shop scheduling problem as a 0–1 integer programming model with the objective of minimizing the total job tardiness along with some assumptions, such as sequence-independent setup and sequence-dependent removal times. They proposed some hybrid geneticbased heuristics to solve the problem in an acceptable computing time. Mosheiov and Oron (2008) addressed batch scheduling problems on an open-shop with m machines and njobs. Identical processing time jobs, machine-independent and sequence-independent setup times are the main assumptions of their problems. The objectives are to minimize the makespan and minimize * Correspondence: [email protected] 2 Department of Industrial Engineering, College of Engineering, University of Tehran, PC: 14399-57131, Tehran, Iran Full list of author information is available at the end of the article © 2012 Noori-Darvish and Tavakkoli-Moghaddam; licensee Springer. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Noori-Darvish and Tavakkoli-Moghaddam Journal of Industrial Engineering International 2012, 8:25 http://www.jiei-tsb.com/content/8/1/25 flow time. They proposed an O(n) time algorithm for the flow time minimization problem. To achieve Pareto-optimal sets for medium to largesized open shop problems using efficient meta-heuristic methods can be necessary and helpful. Naderi et al. (2011) considered an open shop with a set of parallel machines at each stage to minimize the total completion times. They proposed a mixed-integer linear programming (MILP) model for this problem. Moreover, they applied a memetic algorithm to solve the problem. Ahmadizar et al. (2010) addressed a stochastic group shop scheduling problem with known distributions for random release dates and random processing times. They formulated a stochastic programming problem and solved it by the use of an approach being a hybrid of an ant colony optimization (ACO) algorithm and a heuristic algorithm to minimize the expected makespan. Zhang and van de Velde (2010) considered an on-line twomachine open shop scheduling problem with time lags between the completion time and the start time of two consecutive operations of any job. They developed and analyzed the performance of a greedy algorithm to minimize the makespan. Mastrolilli et al. (2010) dealt with a concurrent open shop and proposed a primal– dual 2-approximation algorithm to minimize the total weighted completion times. They also considered several natural linear programming relaxations for the problem. Fei et al. (2010) considered a weekly surgery schedule in an operating theatre. The objectives are to maximize the utilization of the operating rooms, minimize the overtime cost in the operating theatre, and minimize the unexpected idle time between surgical cases. They proposed a column-generation-based heuristic (CGBH) procedure and a hybrid genetic algorithm (HGA) for solving the planning problem as a set-partitioning integer-programming model and the daily scheduling problem as a 2-stage hybrid flow-shop problem, respectively. Matta (2009) developed two original mixed-integer programming (MIP) models (i.e., time-based model and sequence-based model) for the proportionate multiprocessor open shop scheduling problem and proposed a genetic algorithm (GA) to schedule the shop with the objective of the makespan minimization. Seraj and Tavakkoli-Moghaddam (2009) proposed a TS method to solve a new bi-objective mixed-integer mathematical programming model for an OSSP. This model seeks to minimize the mean tardiness and the mean completion time. Panahi and Tavakkoli-Moghaddam (2011) proposed a hybrid method based on multi-objective simulated annealing (SA) and ant colony optimization (ACO) to solve an open shop scheduling problem that minimizes biobjectives, namely makespan and total tardiness. They compared their computational results with a well-known multi-objective genetic algorithm, namely NSGA II. By considering the previous studies, we can conclude that there is only one paper considered the sequencedependent setup time as an assumption for a singleobjective OSSP without any mathematical model and written by Roshanaei et al. (2009). Therefore, in this paper, a bi-objective mixed-integer linear programming (MILP) model is designed with the sequence-dependent setup time as a constraint for the OSSP. To solve medium to large-sized problems, the well-known NSGA-II proposed by Deb et al. (2002) is improved by using a heuristic method in order to achieve good approximate Paretooptimal frontiers. The rest of this paper is organized as follows. The mathematical programming is presented in Section 2. The interactive multi-objective decision making approach is presented in Section 3. Section 4 elaborates the proposed improved NSGA-II. The numerical examples, computational results and performance analysis are given in Section 5. Finally, Section 6 includes the conclusion remark. Mathematical Programming The OSSP considered in this study includes njobs to be processed on at most mmachines. We extend the mathematical model proposed by Low and Yeh (2008) by changing it to a bi-objective model with sequencedependent setup times. Problem assumptions The main assumptions of the presented model are as follows. Each job is to be processed on at most mmachines. The processing sequence of each job is immaterial. A job is not processed by more than one machine simultaneously. Each machine processes at most one operation at a time. All jobs are available at time 0. No preemption is allowed. Machine-dependent and sequence-dependent setup times are considered for each operation. Notations The indices, parameters and decision variables used to formulate the mathematical model are introduced below. Indices NSet of jobs to be processed, N=1;2;...;n fg ;N jj ¼n LSet of machines, L¼1;2;...;m fg ;L jj ¼m. Noori-Darvish and Tavakkoli-Moghaddam Journal of Industrial Engineering International 2012, 8:25 Page 2 of 13 http://www.jiei-tsb.com/content/8/1/25 i,k,g Job indices, (i,k,g= 0,1,...,n),where job 0 is a dummy job. j,h Machine indices, (j,h=1,2,...,m). Parameters MA large positive number. O ij Operation of job ion machine j;8i2N;8j2L. S kij Setup time of job ion machine jimmediately after job k. p ij Processing time of job ion machine j. d i Due date of job i. H j Dependent set up time matrix of machine j. 012 ...n Hj¼ 0 1 2 ... n S01jS02j... S0nj  S12j... S1nj S21j... S2nj ⋮⋮⋮⋱⋮ Sn1jSn2j⋯ 2 6 6 6 6 4 3 7 7 7 7 5 : Decision variables Tsij Starting time of a setup task for operation O ij . T i Tardiness of job i. C max Makespan. Yijh ¼1if 0ij precedes 0ih for jobi 0Otherwise 8i2N;8j;h2L;j6¼ h Xkij ¼1if 0kj precedes 0ij on machine j Otherwise 8k2N[0 fg ;8i2N;i6¼ k;8j2L Zkij ¼1if 0kj precedes 0ij on machine j 0Otherwise 8k2N[0 fg ;8i2N;i6¼ k;8j2L Mathematical model As mentioned at the beginning of Section 2, the model designed in this paper is developed by extending and modifying the mathematical model proposed by Low and Yeh (2008). In this new model, the sequencedependent setup time is used instead of the sequenceindependent setup time considered in the traditional model. The makespan is added to the model as another criterion along with the total tardiness for the minimization purpose. According to these changes and for the validity of the extended model, all of the existing constraints are changed and four new constraints are added to the model. Thus, the bi-objective MILP (BOMILP) model is formulated. It should be noted that S 0ij is the setup time of job ion machine jwhen iis the first job on machine j. Moreover, if job iis the first job on machine j, then Z 0ij = 1. The proposed model is as follows: Min Z1¼X i¼1 Tið1Þ Min Z2¼Cmax ð2Þ s.t. Tsij þX n k¼0;k6¼i ðSkij ZkijÞþ pij≤Cmax ;8i2N;8j2L ð3Þ Tsij þX n k¼0;k6¼0 ðSkij ZkijÞþpij Mð1Yijh≤TsihÞ; 8i2N;8j;h2L;j6¼ hð4Þ Tsih þX n k¼0;k6¼i ðSkih ZkihÞþpih MYijh≤Tsij; 8i2N;8j;h2L;j6¼ hð5Þ Tsij þX n k¼0;k6¼i6¼g ðSkij ZkijÞþpij Mð1XigjÞ≤Tsgj; 8i;g2N;i6¼ g;8j2Lð6Þ Tsgj þX n k¼0;k6¼i6¼g ðSkij ZkijÞþpgj MXigj≤Tsij; 8i;g2N;i6¼ g;8j2Lð7Þ Tsij þX n k¼0;k6¼i ðSkij ZkijÞþpij Ti;8i2N;8j2L ð8Þ Yijh þYihj ¼1;8i2N;8j;h2L;j6¼ hð9Þ Xigj þXgij ¼1;8i;g2N;i6¼ g;8j2Lð10Þ Xkij Zkij≥0;8k2N[0 fg ;8i2N;i6¼ k;8j2L ð11Þ Xkij þZikj≤1;8i;k2N;i6¼ k;8j2Lð12Þ X n i¼1;i6¼k Zkij≤1;8k2N[0 fg ;8j2Lð13Þ X n k¼0;k6¼i Zkij ¼1;8i2N;8j2Lð14Þ Tsij ≥0;8i2N;8j2Lð15Þ Ti≥0;Cmax≥0;8i2Nð16Þ Noori-Darvish and Tavakkoli-Moghaddam Journal of Industrial Engineering International 2012, 8:25 Page 3 of 13 http://www.jiei-tsb.com/content/8/1/25 Xkij;Zkij 20;1 fg ;8k2N[0 fg ;8i2N;i6¼ k;8j2L ð17Þ Yijh 20;1 fg ;8i2N;8j;h2L;j6¼ hð18Þ Two objective functions (i.e., total tardiness and makespan) are shown by Eqs. (1) and (2).Constraint (3) describes the makespan. Constraints (4) and (5) express the relationship between two operations of job ithat do not require being consecutive. The starting time for setup task of operation O ih is greater than or equal to the completion time of operation O ij . Constraints (6) and (7) state the operational sequence of the operations which are processed on the same machine and do not require to be consecutive. The setup task of operation O gj cannot be started until machine jhas finished the processing task of operation O ij . Constraint (8) describes the tardiness for job i. Constraint (9) expresses the order of any two operations of a job; if Y ijh = 1 then Y ihj =0; otherwise, Y ijh =0andY ihj = 1. Constraint (10)expresses the order of any operation pairs (O ij ,O gj ) on the same machine j;ifX igj = 1 then X gij = 0, otherwise, X igj =0 and X gij = 1. Constraints (11) and (12) describe the relationship between X kij and Z kij ;ifX kij = 1, then Z kij =1or 0; otherwise, X kij = 0 and Z kij = 0, by considering the dummy job 0. If X kij = 1 then Z ikj = 0; otherwise, X kij =0 and Z ikj = 1 or 0; in this case, the dummy job 0 is not considered. Because the dummy job cannot be located after any job, it is only used for characterizing the first job on each machine to apply the corresponding relative setup time. Constraint (13) indicates that there is at most one job which can be processed immediately after job kon machine j. if job kis the last job on machine j, then Z kij = 0. Constraint (14) indicates that there is only one job that can be processed immediately before job ion machine jby considering the dummy job 0. Constraint (15) expresses that all jobs should be available for scheduling at time 0.Constraints (16) to (18) define the continuous and binary decision variables, respectively. Interactive MODM Approach To solve the original bi-objective decision making (BODM) problem, an interactive fuzzy programming solution method, called the TH method proposed by Torabi and Hassini (2008), is used. This method is applied to achieve the Pareto-optimal solutions of the presented bi-objective crisp model. Torabi and Hassini (2008) proved that the TH method obtains efficient solutions for the original multi-objective model. According to the characteristics of the given problems, the steps of the TH method are as follows. Algorithm1: TH method Step 1 Calculate the positive ideal solution (PIS) and the negative ideal solution (NIS) for each objective function by solving the corresponding MILP model given below. ZPIS 1¼minX n i¼1 Tis:t:X2FðxÞð19Þ ZNIS 1¼maxX n i¼1 Tis:t:X2FðxÞð20Þ ZPIS 2¼minCmax s:t:X2FðxÞð21Þ ZNIS 2¼maxCmax s:t:X2FðxÞð22Þ where Xis a feasible solution vector consists of all of the continuous and binary variables in the original model and F (x) denotes the feasible area consists of Constraints (3)to(18). Obtainment of the above ideal solutions requires solving four mixed-integer linear programs. In order to reduce the computational time and complexity, we can calculate the PIS for each objective function by solving the corresponding MILP models given in Eqs. (19) and (21). Then, the negative ideal solutions can be determined by the use of the following heuristic rule: ZNIS h¼MaxðZhðx kÞÞ;h¼1;2;k¼1;2ð23Þ It is noted that x hand Zhðx hÞdenote the decision vector associated with the PIS of the h-th objective function and the corresponding value of the h-th objective function, respectively (Torabi and Hassini, 2008). The related results are shown in Table 1. Step 2 For each objective function, determine a linear membership function by: μZ1ðXÞ¼ 1;Z1<ZPIS 1 ZNIS 1Z1 ZNis 1ZPIS 1 ;ZPIS 1≤Z1≤ZNIS 1 0;Z1>ZNIS 1 8 > > < > > : ð24Þ μZ2ðXÞ¼ 1;Z2<ZPIS 2 ZNIS 2Z2 ZNis 2ZPIS 2 ;ZPIS 2≤Z2≤ZNIS 2 0;Z2>ZNIS 2 8 > > < > > : ð25Þ Figure 1 illustrates the graph of these membership functions. Noori-Darvish and Tavakkoli-Moghaddam Journal of Industrial Engineering International 2012, 8:25 Page 4 of 13 http://www.jiei-tsb.com/content/8/1/25 Step 3 Transform the BOMILP model into an equivalent single-objective MILP using the following auxiliary formulation. Max λðXÞ¼γλ0þð1γÞXhθhμZhðXÞð26Þ s.t. λ0≤μZhðXÞ;h¼1;2ð27Þ X2FðXÞð28Þ λ0;γ2½0;1ð29Þ According to the two objective functions considered in the problem, Constraint (27) can be written by: ZNIS 1X n i¼1 Ti≥λ0ðZNIS 1ZPIS 1Þð30Þ ZNIS 2Cmax ≥λ0ðZNIS 2ZPIS 2Þð31Þ where μZhðXÞis the satisfaction degree of the h-th objective function and λ0¼minhμZhðXÞ  is the minimum satisfaction degree of objectives. Also, θ h denotes the importance level of the h-th objective function such that Phθh¼1;θh>0 . The θ h parameters are determined linguistically by the decision maker based on her preference. Moreover, γis the coefficient of compensation. By changing the values of this parameter in the interval [0,1], the TH method can obtain both unbalanced and balanced compromised solutions. It means that for higher values of γ, the solution method results bigger lower bound for the satisfaction degrees of objectives (λ 0 ) for a given sample example. These solutions are balanced compromised solutions. These kinds of solution can be more appropriate when the importance levels of all objective functions are equal. On the other hand, for lower values of γ, the solution method resulting solutions with bigger satisfaction degrees for some objectives with higher importance levels than others. These solutions are unbalanced compromised solutions. These kinds of solutions can be more appropriate when the importance levels of the objective functions are different. Step 4: Solve equivalent single-objective MILP model using the given coefficients (θ h ,γ). If the decision maker is satisfied with obtained efficient compromised solution, then stop; otherwise, change the value of some controllable parameters and then go to Step 2. A set of non-dominated solutions as the Paretooptimal solutions can be obtained for the BOMILP model by changing the values of controllable parameters of the MILP model given in Eqs. (26) to (31). As mentioned at the beginning of Section 3, it is proved that the TH method achieves efficient solutions for multiobjective optimization problems. It means that each optimal solution of the single-objective model given in Eqs. (26) to (31) is a Pareto-optimal solution of the original BOMILP model. Proposed Method The OSSP studied in this paper is the NP-hard problem. Thus, to solve medium to large-sized problems considered in this paper, an improved non-dominated sorting genetic algorithm II (NSGA-II) is proposed. Noori-Darvish and Tavakkoli-Moghaddam (2011) used the original NSGA-II for an OSSP. NSGA-II belongs to a class of multi-objective evolutionary algorithms (MOEAs). In most problems, it is able to find much better spread of solutions and better convergence near the true Pareto-optimal front compared to two other elitist MOEAs, namely Pareto-archived evolution strategy (PAES) and strength-Pareto evolutionary algorithm (SPEA), which pay special attention to creating a diverse Pareto-optimal front (Deb et al., 2002). Solution Representation In our proposed method, a chromosome is an operation-based array or permutation list that has been a typical solution representation studied in OSSPs. As illustrated in Figure 2, the permutation list is a single-row array consisting of n×moperations. By using this encoding scheme, the basic assumptions of OSSPs are satisfied, in which the processing sequence of each job is immaterial and a job is not processed more than once by one machine. In this representation, operations are Table 1 Payoff results Z1Z2 χ 1ZPIS 1ZNIS 2 χ 2ZNIS 1ZPIS 2 1 Figure 1 Linear membership function for Z 1 (Z 2 ). Noori-Darvish and Tavakkoli-Moghaddam Journal of Industrial Engineering International 2012, 8:25 Page 5 of 13 http://www.jiei-tsb.com/content/8/1/25 listed in the relative order by which they are scheduled. pdenotes the position of an operation in the list. To decode the list presented in Figure 2 and obtain its equivalent solution in the solution space, we start from the beginning of the list and schedule the first operation. Then, according to the corresponding order of operations in the list, we schedule other operations one by one. It means, job 5 is first scheduled on machine 3. Then, job 2 is scheduled on machine 4 and so on. According to the rule that each operation is scheduled at the earliest time, there are two variables, namely 1) Max {C ij } for all ias the maximum completion time of jobs on machine jand 2) Max {C ij } for all jas the maximum completion time of job i. The starting time and completion time of each operation are computed by: Tsij ¼Maxf0;Max 8ifCijg;Max 8jfCijgg ð32Þ Cij ¼Tsij þX n k¼0;k6¼i ðSkij ZkijÞþpij ð33Þ Initialization A good initial population can help an evolutionary algorithm to start multi-objective optimization from good solutions in the solution space and have a better performance that will be discussed in the next section. In our algorithm, a heuristic method is applied to generate the initial set of solutions or population. It first constructs a set of random sequences of operations as chromosomes (as many as the population size). By using the swapping of two adjacent or non-adjacent operations (see Figures 3 and 4), it generates all the feasible sequences in the neighborhood of each permutation list. For each chromosome, the best feasible solution obtained by these neighborhood searches is a member of the initial population. According to the dominance concept, the best solution is an efficient solution of the Pareto-optimal frontier. The set of Pareto-optimal solutions dominate any other solutions in the feasible area and improve all the objective functions simultaneously. It should be noted that k,k' = 1,2,...,n×m, as shown in Figures 3 and 4. Crossover A procedure proposed by Low and Yeh (2008) is used for the crossover operator. Algorithm 2: Crossover operator Step 1 Select two cut points randomly along the positions of the strings for each pair of parent strings (A and B). So, each permutation list is divided into three sections, called substring 1, substring 2 and substring 3. Step 2 Exchange substring 1 of parent string A and substring 3 of parent string B. Similarly, exchange substring 3 of parent string A and substring 1 of parent string B. Do not change the elements of substring 2 of each parent string. Therefore, two proto-children are generated which are named A 0 and B 0 . Step 3 Legalize two generated offsprings A 0 and B 0 .In order to do this, remove the elements of substring 1 and substring 3 of offspring string A 0 , which are the same as substring 2, and then replace them with corresponding elements of substring 2 of offspring string B 0 . Also, do the same procedure for offspring string B 0 . Thus, two offspring strings (i.e., A 0 and B 0 ) are produced. Mutation As depicted in Figure 5, the insertion operator is applied for the mutation operator. In the permutation shown in this figure, k,k’= 1,2,....,n×m. p=n×m…p=4p=3p=2p=1 Oij …O11 O34 O24 O53 Figure 2 Operation-based representation/permutation list. p=n×m ………. p=k+1 p=k ……… Oij ……... O” ij O’ij …….. Figure 3 Swap of two adjacent operations in a permutation. p=n×m ……….. p=k’ ……….. p=k ………. Oij ………. O” ij ………. O’ij …….. Figure 4 Swap of two non-adjacent operations. p=n×m ……….. p=k’ ……….. p=k ………. Oij ………. O” ij ………. O’ij …….. Shift Figure 5 Insertion operator for mutation. Noori-Darvish and Tavakkoli-Moghaddam Journal of Industrial Engineering International 2012, 8:25 Page 6 of 13 http://www.jiei-tsb.com/content/8/1/25 Elitism and selection operator Deb et al. (2002) designed an elitist NSGA including “Fast non-dominated sorting approach”and “Diversity Preservation”for solving multi-objective problems. In order to rank solutions (individuals) and sort them into different non-dominated frontiers, an approach called “Fast non-dominated sorting”is applied. By the use of this approach, a domination count (n p ) is calculated for each solution to specify its non-dominated frontier. The domination count of all solutions in the first nondominated frontier (i rank =1) is equal to zero. At the end of the multi-objective optimization process, solutions with n p =0, dominate all other solutions in the solution space. It means they are the best global solutions. This approach results better convergence near the true Pareto-optimal frontier. Along with the convergence to the Pareto-optimal set, it is also desired that an EA maintains a good spread of solutions in the obtained set of solutions (Deb et al., 2002). An operator called crowded-comparison operator is proposed for “Diversity Preservation”. By considering the value of crowding distance of each solution (i distance ) in its frontier, this operator guides the selection process at the various stages of the algorithm toward a uniformly spread-out Pareto-optimal frontier. When the improved NSGA-II is iterated as many as the pre-specified iterations (i.e., maximum number of generations), the multi-objective optimization process is terminated (Noori-Darvish and Tavakkoli-Moghaddam, 2011). Results and discussion Several numerical examples in small to large sizes are generated randomly by the use of a classic approach of the literature in order to examine and analyze the validity and efficiency of the mathematical model presented in Section 2, the TH method presented in Section 3, and the performance of the improved NSGA-II presented in Section 4. The small-sized problems are solved exactly by the use of the Lingo 9 software and the results of the TH method are analyzed. If the number of jobs and the number of machines are more than 4, the sizes of the problems are medium to large. In these conditions, even after several hours of running, Lingo cannot yield an optimal solution in a reasonable time. Thus, the improved NSGA-II is applied to solve these kinds of problems. Finally, the efficiency of this algorithm is first compared with the efficiency of a well-known multi-objective genetic algorithm, namely SPEA-II, by using the design of experiments (DOE). Then, the proposed algorithm is compared with the traditional NSGA-II. These algorithms are coded in Turbo C++4.5 on a PC (Main board P4, CPU E5200 2.5GHz, 6M Cache, RAM 2GB BUS 800). SPEA-II The strength Pareto evolutionary algorithm II (SPEA-II) proposed by Zitzler, et al. (2001) is specially designed for multi-objective optimization problems. This algorithm includes some special features (i.e., fitness assignment strategy) considering a number of solutions that dominate each individual and a number of solutions that each individual dominates the “nearest neighbor density estimation technique”that yields a density value for each solution, and the “archive truncation method”that preserves a specified number of solutions in the external non-dominated archive. In our SPEAII, the solution representation is the permutation list and the initial population is generated randomly. Moreover, the binary tournament selection procedure is applied. When the SPEA-II is iterated as many as the pre-specified number of iterations (i.e., maximum number of generations), the multi-objective optimization process is terminated. Generating numerical examples In this section, we use a classical approach in the literature proposed by Loukil et al. (2005) to generate several numerical examples of small to large sizes Table 2 Parameters of the NSGA-II Population size Maximum number of generations Crossover rate Mutation rate 100 300 0.95 0.02 Table 3 Parameters of the SPEA-II Population size Maximum number of generations Selection rate Crossover rate Mutation rate 100 400 0.8 0.8 0.2 Table 4 Values of parameters for numerical instances Numerical instances RT a 1 0.8 b 0.6 0.6 c 0.2 0.4 Table 5 Values of positive and negative ideal solutions Z 2 Z 1 Sample NIS PIS NIS PIS example 119 95 193 130 4×3-a 103 99 114 107 4×3-b 112 104 52 26 4×3-c 141 119 383 283 4×4-a 151 146 310 196 4×4-b 145 135 249 136 4×4-c Noori-Darvish and Tavakkoli-Moghaddam Journal of Industrial Engineering International 2012, 8:25 Page 7 of 13 http://www.jiei-tsb.com/content/8/1/25 randomly. The processing times and due dates are uniformly distributed in the intervals [0,100] and P1TR 2  ;P1TþR 2  , respectively. where P¼ðmþn1Þ Pand  P¼Pn i¼1Pm j¼1pij=ðnmÞ which is the mean values of the processing times. Parameters Rand Ttake their values in the sets {0.2,0.6,1}, {0.4,0.6,0.8},respectively. Also, the setup times are random variables between 10 and 50. Parameter setting Various sets of controllable parameters and different sizes of problems are considered. Then, many experiments are designed and run by using them. The effective values in terms of solution quality are determined. The parameter of the TH method are as follows. Parameter γ takes its values in the set {0,0.1,...,1}. Also, it is assumed that the preference information corresponding to the importance levels of the objective functions are specified linguistically by the decision maker as: θ 1 =θ 2 and θ 1 > θ 2 . So, the values of these parameters in the first condition are θ 1 =θ 2 =0.5, and in the second condition are θ 1 =0.8 and θ 2 =0.2. It should be mentioned that the values of controllable parameters (i.e., Rand T) of the generating instances approach are given in the following subsections. Tables 2 and 3 depict the values are set for the parameters of the NSGA-II and SPEA-II, respectively. Moreover, each example is solved 10 times independently. Solving small-sized problems Two kinds of sample examples in small sizes with 4 jobs and 3 machines (i.e., 4×3), and 4 jobs and 4 machines (i.e., 4×4) are considered. For each type of these examples, three numerical instances are generated randomly. Tables 4 and 5 represent the values of controllable parameters and the values of positive ideal solutions (ZPIS i) and negative ideal solutions ( ZNIS i) for each numerical instances, respectively. Using the obtained values of the PIS and the NIS of each objective function, the final MILP model is exactly solved for each numerical instances by the Lingo 9 software. Table 6 illustrates the results of example 4×4-b. The obtained resuts are analyzed as follows. In more cases with different importance levels of objective functions, the TH method performs correctly and efficiently. It means according to the importance levels , the solutions found by this method are unbalanced compromised solutions. However, in some cases, the TH method does not perform well and the satisfaction degree of the objective function with lower importance is more than that of the objective function with higher importance. When the importance levels are equal, the satisfaction degrees of objective functions are very close to each other, in some cases. Thus, the method approximately finds balanced compromised solutions. According to the discussion in Section 3, each optimal solution of the final MILP model is an efficient solution to the BOMILP model. Thus, in all cases, by changing the values of controllable parameters (i.e. γand θ) two Pareto-optimal solutions are found by the TH method. Figure 6 indicates the Pareto-optimal frontier of example 4×4-b. Table 6 Computational results of example 4×4-b θ 1 = 0.8 , θ 2 = 0.2 θ 1 = 0.5 , θ 2 = 0.5 γ μZ2μZ1Z 2 Z 1 μZ2μZ1Z 2 Z 1 0.4000 0.9298 149 204 0.4000 0.9298 149 204 0 0.4000 0.9298 149 204 0.6000 0.5614 148 246 0.1 0.4000 0.9298 149 204 0.6000 0.5614 148 246 0.2 0.4000 0.9298 149 204 0.6000 0.5614 148 246 0.3 0.4000 0.9298 149 204 0.6000 0.5614 148 246 0.4 0.6000 0.5614 148 246 0.6000 0.5614 148 246 0.5 0.6000 0.5614 148 246 0.6000 0.5614 148 246 0.6 0.6000 0.5614 148 246 0.6000 0.5614 148 246 0.7 0.6000 0.5614 148 246 0.6000 0.5614 148 246 0.8 0.6000 0.5614 148 246 0.6000 0.5614 148 246 0.9 0.6000 0.5614 148 246 0.6000 0.5614 148 246 1 T= 0.6 R= 0.6 145 146 147 148 149 150 151 152 0 50 100 150 200 250 300 350 Z 2 Z1 Ideal point 1 Ideal point 2 Test problem Figure 6 Pareto-optimal frontier of example 4×4-b. Table 7 Characteristic of medium to large-sized examples Number of machines Number of jobs Representation Sample example 5 10 10×5 1 7 14 14×7 2 9 20 20×9 3 Noori-Darvish and Tavakkoli-Moghaddam Journal of Industrial Engineering International 2012, 8:25 Page 8 of 13 http://www.jiei-tsb.com/content/8/1/25