CMS scheduling problem considering material handling and routing flexibility

Cell manufacturing as an application of group technology increases the flexibility and efficiency of the production. Cell scheduling problem, one of the subjects in cell manufacturing, has not been widely studied by researchers compared with other problems in cell manufacturing. In spite of great importance of material handling in cell scheduling, it has not been paid enough attention by researches. In this paper, a new mathematical model for cell scheduling problem considering material handling time and routing flexibility is proposed. The proposed model belongs to the mixed-integer nonlinear programs (MINLP). A linearization procedure is proposed to convert the MINLP to an integer program (IP) in order to develop more powerful optimization tools. Furthermore, a simulated annealing-based heuristic is developed to solve the large-size problems.


Introduction
In the recent decades, production systems have changed over due to the increase of competition in markets. In the past, production volume and finished costs of goods were the two major determinant competitive factors; thus, the tendency to use flowshop production systems was pervasive among the companies. Nowadays, other factors such as variety in products and swift response to the market's demand are of higher importance. Consequently, production systems with higher level of flexibility have been developed and applied by both practitioners and academics. Job shop can be considered as a clear-cut example of such systems.
The aforementioned changes along with the development of new businesses have created a novel environment for manufacturing and competing in markets. In such an environment, reducing production costs and make-span, increasing flexibility, and precipitating reactions to the market's needs are turning into the critical competitive advantages of companies. In the 1970s, through highly competitive atmosphere among the companies in the US, a few modern management concepts such as just in time (JIT) and group technology (GT) came into existence. Cell manufacturing system (CMS) is considered as an application of GT, in which the whole machines are divided into a number of distinctive groups according to the similarities in processing the assigned parts. The groups are called cells. The set of parts processed in each cell is called a part family. Cell formation, cell layout, production planning, and cell scheduling are the major issues in CMS (Solimanpur et al. [25]). Unlike the cell formation problems, cell scheduling has not been widely studied by researchers. Moreover, in spite of the great importance of material handling issue, it has not been considered in the most of cell scheduling problems.
As Nomden and Van der Zee [19], routing flexibility provides the possibility to choose from among a number of machines to execute an operation. We consider routing flexibility in terms of alternative machines available for a product family. On the other hand, there will be alternative production routes per product family. This issue can increase the flexibility of a manufacturing system; however, it can increase the complexity of the problem.
In the current research, a new mathematical model for a cell scheduling problem is proposed. The main contributions of the proposed model are explained as follows: -There are some papers that have considered routing flexibility in cellular manufacturing problem such as Kioon [14]. To the best of our knowledge, this issue has not been directly integrated with CMS scheduling problem while it may have great influence on the sequencing of the jobs. -Material handling time is also one of the subjects, which is undermined in the literature of CMS scheduling. In the proposed model, the time is determined based on the characteristics of the part and position of the machines. -A number of the researches in this sphere assume that all the operations of a part family on a particular machine have to be done consecutively and without any interruption [16,25]. This assumption leads to minimizing the total setup times and may be interpreted as a proper policy for minimizing the makespan. This assumption reduces the solution space and as a result can decrease the computational time. However, it is clear that in many cases, consecutively, performing operations of a part family is not optimal and may impede the solving algorithm from finding the global optimum. In the given model, a machine can process the parts without any restriction on the sequence considering the setup times. This approach extends the applicability of the proposed model to the job shop environment scheduling.
The complexity of scheduling problems and the nonlinearity of the models hinder the applicability of the regular optimization tools. Therefore, we alleviate this problem by converting the nonlinear model into a linear one in order to implement more efficient optimization tools. Finally, a simulated annealing (SA) based heuristic is proposed in order to tackle the large-size instances of the problem. This paper is organized as follows: a literature review of CMS scheduling is presented in Section 2. In Section 3, the problem description and the mathematical model are given. The linearization process is demonstrated in Section 4, and we introduce the proposed SA-based heuristic in Section 5. Section 6 gives the computational results. In Section 7, conclusions and further research ideas are given.

Literature review
Generally, a significant portion of the literature in CMS is devoted to the flowline-based cellular manufacturing. In this sphere, it is assumed that all the parts follow a similar sequence of processes; however, there may be some missing operations on some machines. Since the proposed model in this paper can be considered as a submodel in cell scheduling, we directly go through the corresponding literature.
There are a few researches focusing on the minimization of tardiness of jobs in cell scheduling problem. Parthasarathy and Rajendran [20] study the problem of scheduling in flowshop and flowline-based manufacturing cell (or simply, a cell) with the objective of minimizing mean tardiness of jobs. A heuristic algorithm, based on the SA technique, is developed. The proposed SA algorithm with two schemes is evaluated against the existing heuristics that seek to minimize mean tardiness of jobs. The results of the computational evaluation reveal that the SA algorithm with two schemes performs better than the existing rival heuristics. Rajendran and Ziegler [22] study the addressed problem with the objective of minimizing the sum of weighted flowtime and weighted tardiness of jobs. Initially, heuristic preference relations are developed considering lower bounds on the completion times, operation due dates, and weights for holding and tardiness of jobs; then, a heuristic algorithm is proposed making use of the heuristic preference relations. Two more heuristic algorithms are developed, implementing an improvement scheme on the solution given by the first heuristic algorithm.
There are a number of successive researches on the scheduling of part families and jobs within each family in a flowline or flowshop (all jobs available at time zero, different job availability times known a priori) manufacturing cell with sequence-dependent family setups times. Schaller et al. [23] study the addressed problem. The objective is to minimize the makespan while processing parts (jobs) in each family together. França et al. [6] study the problem in a flowshop manufacturing cell. Two evolutionary algorithms including a genetic and a memetic algorithm with local search are proposed and empirically evaluated as to their effectiveness in finding optimal permutation schedules. Hendizadeh et al. [9] present various tabu search (TS)-based meta-heuristics for the problem to minimize makespan. Concepts of elitism and the acceptance of worse moves inspired from simulated annealing are considered in the proposed meta-heuristics to improve intensification and diversification. The effectiveness and efficiency of the proposed heuristics are compared against the best rival meta-heuristic and heuristic algorithms. Lin et al. [16] declare that almost all published studies in this regard focus on using permutation schedules to deal with flowline manufacturing cell scheduling problem with sequencedependent family setup times. To explore the potential effectiveness of treating this argument using nonpermutation schedules, three prominent types of meta-heuristics including SA, genetic algorithm (GA), and TS are proposed and empirically evaluated. The experimental results demonstrate that the improvement made by nonpermutation schedules over permutation schedules for the due date-based performance criteria were significantly better than that for the completiontime-based criteria. Bouabda et al. [2] address the permutation flowline manufacturing cell with sequence-dependent family setup time problem with the objective to minimize the makespan criterion. A cooperative approach including a GA and a branch and bound (B&B) procedure is developed. The latter is probabilistically integrated within the GA.
There are a number of other dispersed researches related to the cell scheduling problem which are reported here. Sridhar and Rajendran [27] develop a multiobjective model minimizing the makespan, the total flow time, and the machine idle time. They utilize a GA-based heuristic to solve the proposed problem. Solimanpur et al. [25] declare that the scheduling problem in a cellular manufacturing environment is treated as group scheduling problem, implying that all parts in a part family to be processed in the same cell. However, there could be some exceptional parts, which need to visit machines in the other cells. This fact limits the applicability of group scheduling approaches. The addressed research study the scheduling of manufacturing cells in which parts may need to visit different cells. A two-stage heuristic is proposed to solve this problem. These stages are termed as intracell scheduling and intercell scheduling. Through intracell scheduling, the sequence of parts within manufacturing cells is determined. In intercell scheduling, however, the sequence of cells is obtained. Sidhartha and Canel [4] focus on the problem of scheduling batches of parts in a flexible manufacturing system (FMS). Due to the use of serial access material-handling systems in many FMSs, the problem is modeled for a multicell FMS with flowshop characteristics. A B&B solution method is developed in order to solve the problem. Das and Canel [4] propose a new multicell scheduling model. They also developed a B&B method to solve the problem. Although they developed a model assuming multiple cells, the processes of parts on machines were not distinguished, and only the process times of parts were given on each cell. Considering this assumption, they dominated the characteristics of flowshop on their proposed model. Venkataramanaiah [30] develop a SAbased heuristic to minimize makespan, flow time, and machine idle time in cellular manufacturing environment.
Kesen and Gungor [13] design a GA-based heuristic approach for job scheduling in virtual manufacturing cells (VMCs). In a VMC, machines are dedicated to a part as in a regular cell, but machines are not physically relocated in a contiguous area. Cell configurations are therefore temporary, and assignments are made to optimize the scheduling objective under changing demand conditions. This research considers the case where there are multiple jobs with different processing routes. There are multiple machine types with several identical machines in each type. Weighted makespan and total traveling distance are considered as the scheduling objectives. Gholipur-Kanani et al. [8] develop a new model which minimizes the number of intercellular movements as well as intracellular movements, makespan tardiness, and sequence-dependent set-up time, however, not considering routing flexibility and material handling time. In addition, they posed the characteristics of the flowline-based cellular manufacturing to their model assuming that the process of each part family cannot be interrupted by other parts. This assumption can increase the machines' idle time and consequently the makespan significantly. Chalapathi Pasupuleti [3] declares that once the cellular manufacturing system is designed, scheduling of jobs is essential for the day-to-day production in the machine cells. A methodology for prioritizing the parts, as well as preparing the total schedules in cellular manufacturing system, is proposed in this research. It takes into account the processing sequences of the jobs, processing and setup times, and due dates. The method presents the sequence of parts to process on each machine and the total schedules for all the operations of the parts. Detailed reviews of flowline-based cellular manufacturing scheduling problem can be found in Allahverdi et al. [1].
Although there are a number of researches in cell scheduling, but there is not any research considering material handling time and routing flexibility. Jensen et al. [10] believes that for a parallel shop, full routing flexibility can be unfavorable because this results in extensive set-up efforts; instead, the tradeoffs between routing flexibility and setup efficiency must be made carefully. Sheikhzadeh et al. [24] declare that routing flexibility is useful in the cases where either the number of machines or length of setup times is high. Tsubone and Horikawa [29] were of the opinion that increasing routing flexibility always makes progress in the performance, despite the presence of family set-ups. Garavelli [7] believes that this holds especially for high shop loads. Diallo et al. [5] expresses that routing flexibility allows the machines usage factor and the throughput time to be improved. The aforementioned authors seem to unite in that little investment in routing flexibility results in significant performance improvements.
Since we have developed a few meta-heuristic-based heuristics in this research, it is good to point that the application of meta-heuristic algorithms to the design and manufacturing problems is growing in the recent years. Yildiz [31], Yildiz [32], and Yildiz [33] are just three sample researched in this regard. Since the current problem in this research deals with routing and material handling, vehicle design and scheduling can be added to the problem. The addressed heuristics are also applicable to the vehicle design problems using the research by Yildiz and Solanki [34]. According to the published researches, scheduling is an NP-hard problem. As a result, researchers have tried to develop heuristic and meta-heuristic algorithms for cell scheduling problems (such as [11,17,[25][26][27][28]). The reader can refer to Ponnambalam [21] for an exhaustive review.

Problem description
A cellular manufacturing problem including n parts (i= 1,…,n) and m machines ( j=1,…,m) is considered. A number of predetermined processes should be done on each part. Notation p represents the processes ( p=1,2,…,P i ) and P i represents the number of processes to be done on part i. Operations required for manufacturing parts are determined in the production planning. Based on the similarities between production processes, parts are divided into K part families (k= 1,…,K). Each part family requires a common setup time on each machine. In addition to a specific setup time for each part family, each part may need a separate setup time which is assumed to be considered in the corresponding process time. The objective function is to minimize the total processing times of all parts.

Assumptions
& A common setup time is required for part families on each machine. & Setup time for each individual part, when needed, is considered as a part of the corresponding process time. & Some processes for a particular part can be done on more than one machine; this is the concept of the routing flexibility. & It is necessary to process all parts of a part family on a particular machine without any interruption from other part families.

Decision variables
The parameters A ij , H i ,T ij , PC ik , D, S jk are inputs. In the following, extra explanations are given to clarify some of the aforementioned parameters: & A ij is an n×m matrix which assigns processes of parts to the related machines. In the proposed model, it is a nonzero-one matrix. For instance, assume a problem with three parts and four machines. A ij can be as the following matrix: Obviously, the matrix illustrates that part 1 needs two processes; the first process can be done only on machine 1, while the second one can be done on both machines 2 and 3. Other parts can be interpreted the same way.
& An example of the matrix PC ik is illustrated in the following: In this example, three parts are divided into two part families. Parts 1 and 2 form one family while part 3 forms the other family. Apparently, it is a zero-one matrix with only one in each line.
The objective function in (1) represents the makespan given by constraint (9) in which G represents the maximum of the completion times of all processes of all parts. Constraint (2) ensures that the process p of part i is done on exactly one machine. As constraint (3), at most one process, of part i can be done on machine j.
Constraints (4), (5), and (6) provide lower bounds for makespans. The first process of part i is the trigger of its total makespan. From the second process onwards, the makespan of the former process of part i, traveling time between machines, and the waiting time for the idle machine affect the total makespan. Consequently, constraint (4), which is just written for the first process of part i, only considers the process and setup times. The first term represents the sum of process time of part I; the second term is the setup time added to the first term if parts i and i′, the previous part on machine j, do not belong to the same part family. The third term is added as a setup time if part i is the first part of its family which is processed on machine j.
Constraint (5), which is written for all processes of the parts except for process 1, takes into account the traveling time between machines by term four and makespan of the former process of part i by term 5 as well as the explained three terms in constraint (4) for the first process of each part. If part i reaches machine j while it is in use by another part, it should wait. This issue is considered in the last term of constraint (6). In fact, the makespans of part i must be greater than the maximum of the values given by constraints (5) and (6) from the second process onwards. If part i reaches machine j while it is in use, the value in constraint (6) is greater than that of constraints (5), while, when machine j is free and is immediately used by part i, the value in constraint (5) is greater than that of constraint (6). Constraint (7) indicates that part i ' can be the former part of part i on machine j, if and only if each of part i and part i ' has a process on machine j. Constraint (8) satisfies the routing flexibility assumption and restricts the machines that can be used for each process based on the PMIM matrix. Constraint (10) illustrates the simple fact that there is no former part for the first part on machine j. The fact that a particular machine cannot operate on more than one part simultaneously is guaranteed by constraints (11) and (12). Constraint (11) guarantees at most one part can be processed after part i on the same machine, while constraint (12) ensures that not more than one part can be operated just before part i.

Numerical example
A clear-cut numerical example is presented in this section in order to clarify the performance of the proposed model. Imagine a CMS problem with four parts, four machines, and with two cells. Input parameters are given as follows: We utilize Lingo solver to solve this example using B&B method. The results are given in Table 1.
The first part family consists of parts 1 and 2, while parts 3 and 4 form the second part family. Figure 1 depicts the final solution schematically. It indicates operational sequences with time consideration. The left side number in the container represents the number of the part, and the right side number represents the number of process done on the part. In addition, movements are drawn using the spiral arrows for the parts. As it is clear from the figure, the first process of parts 2 and 4 on machines 1 and 3 are the starter of the whole process. It takes two units for set-up time and five more units for the process time, so the first processes of these two parts are completed at the seventh minute. Part 2 is then moved to the second machine; considering the handling, setup, and processing time, the second operation of this part is completed at 19th minute. Totally, 24 min is required to complete all the processes of all the parts.

Optimization of nonlinear problem
As it was mentioned in Section 3.1, the presented nonlinear model is solved using Lingo solver on a PC with 2 GHz processor and 2 GB RAM for different sizes of problems. Table 2 gives the size of the problems and the corresponding CPU runtime.
In Table 2, NP, NM, and NC represent the number of parts, number of machine, and the number of part families, respectively. Since the solver uses B&B method, F best represents the best feasible objective function value while F bound represents the lower bound of the objective function value. The final solution is optimum if F bound equals to F best . We stop the solver when the computation time reaches 14,400 s. The solver is capable to find the optimal solution for the first two problems. As the size of the problem increases, the solver presents a feasible solution and a lower bound for the problem in a reasonable processing time (the third problem). In the fourth problem, the solver is not even capable of finding a feasible solution within 4 h, although this problem is considered as a small-size problem.

Linearization of the model
The unsatisfactory performance of the aforementioned solver regarding bigger sizes of the problem motivated us to modify Table 1 Final solution for the given numerical example the structure of the problem in such a way as to implement more efficient optimization tools. Some terms in the constraints are nonlinear which cause intricacy in the model. In this section, the linearization of the mathematical model using some new variables is presented.
ÞÂS jk can be linearized by the following terms: st : Proof Consider the following statements: In the first statement, V ii 0 jp ≥ 1 . In as much as this is a minimizing problem and PC ik Â 1−PC i 0 k À Á Â S jk ≥ 0 , V ii 0 jp will have the least possible value which is 1; therefore nothing will change in the optimal solution. In the other cases, V ii 0 jp ≥ 0 and with the same logic V ii 0 jp would be equal to zero and the optimal solution will not change.
The same explanation can be given for exchanging the following preposition for linear terms: can also be written as follows: The second nonlinear term can be linearized as follows: st : Proof The following cases are possible: In the first statement, B ii 0 jp ≤ 1 . Since the second nonlinear term has a negative sign and the problem is of minimizing, thus B ii 0 jp will have the maximum possible quantity which is 1. For the other cases, either B ii 0 jp ≤ 1 2 or B ii 0 jp ≤ 0 , B ii 0 jp can only have zero value because it is a binary variable.
linearized by the following terms: st : Proof As we have hypothesized in preposition1: Owing to this, the nonlinear term in this section will turn into the following terms: X j In the next step, V ii 0 jp Â x i 0 jp 0 is replaced with U ii 0 jpp 0 along with the constraints below: st : Consider the following six statements below: As the first statement, V ii 0 jp ¼ 1 and x i 0 jp 0 ¼ 1 . As a result, U ii 0 jpp 0 ≥ C i 0 p 0 . Due to the minimizing nature of the model, U ii 0 jpp 0 will have the least value, which is C i 0 p 0 .
In the other cases, U ii 0 jpp 0 ≥−M . M stands for a very big positive number. Since the model is of minimizing and U ii 0 jpp 0 ≥ 0 , U ii 0 jpp 0 will get zero value.
earized by the following terms: Proof According to (37), z ii 0 j can only be equal to 1 when both x ijp and x i 0 jp 0 are equal to 1. On the other hand, according to the rest of the constraints, z ii 0 j would never be equal to 0 when both x ijp and x i 0 jp 0 are equal to 1.
Proposition 6: G=Max{C ip } can be linearized by replacing it with the following set of constraints: Proof According to (38), G is greater than all the value of C ip , for every i and p. Since the problem has a minimization form, G is equal to the maximal value of completion times.

Optimization of the linear model
Linearization of the model brings up the possibility of implementing more specialized tools for integer linear programming. In this section, we solve a set of problems with CPLEX and the results are demonstrated in Table 3. The same PC as expressed in Section 3.2 is used to solve the problems. As given by Table 3, we have been able to find the optimum solution for a problem with 25 machines and 35 parts which are divided into 6 part families (problem no 10). Regarding the optimization results of the nonlinear model, this is a huge improvement in term of CPU run time. For the 11th problem, the CPLEX was unable to find the optimum solution within 14,400 s.

Solution algorithms
Clearly, the considered problem is NP-hard. In Section 4, we tried to modify the problem in such a way as to find the optimal solution for bigger sizes of the problem. However, due to the NP-hard nature of the problem, solving the big sizes of the problem is still impossible within a reasonable time.
Therefore, heuristic algorithms are applied for this situation. Venkataramanaiah [30] indicated the promising performance of SA in CMS scheduling problems. Lin et al. [16] indicated the superiority of SA over GA and TS in this context. These researches motivated us to select SA algorithm from among all the existing options in order to solve the under study CMS scheduling problem.

Simulated annealing
SA is an imitation of the annealing process, in which a material is heated and then allowed to cool very slowly until it reaches its most regular possible crystalline state with corresponding minimum energy [12]. Metropolis et al. [18] was the first research, which suggested the structure of this algorithm. In optimization, SA was used by Kirkpatrick et al. [15] for the first time [12]. The key feature of SA is that it provides a condition by which we can escape local optima by allowing hill climbing moves (i.e., moves which worsen the objective function value) in the hope of finding a global optimum [12]. Kirkpatrick et al. [15] indicated that the Metropolis algorithm could be applied to discrete optimization problems by defining feasible solutions as states, and using an objective function to represent the change in energy between states.

Implementation of SA
In this section, the utilization of SA in order to solve the current problem is explained. Figure 2 depicts the different stages of the designed SA-based heuristic.
Different steps are explained as follows:

Initial solution generation
At the first step, the structure of the solution should be defined.
As it is depicted in Fig. 3, an n×m matrix is proposed to represent each solution of the problem.  Fig. 2 Stages of the designed SA-based heuristic Fig. 3 Initial solution structure In this structure, l sj represents the number of the part whose process is supposed to be done on machine j with order s. For instance, if l 23 is equal to 4, it means that the second part whose process is to be done on machine 3 is part 4. Two methods are proposed in order to generate the solution as follows: & Random generation: It is a simple method in which the initial solution is generated randomly. & Using shortest processing time (SPT) method: Previous researches reveal that using a proper initial solution can improve the performance of a meta-heuristic algorithm; consequently, SPT method was applied to generate a convenient initial solution for the SA-based heuristic. According to SPT, the total sum of processing times of each part family is calculated initially; then, part families are ranked in ascending order based on their processing times. Part families' ranking determines the sequence of processing. The part family with the shortest process time will be of higher priority. The sequence of the parts in each part family is also determined in the similar way, in which the sum of processing time of each part is computed and they are ranked in an ascending order.

Modification of the solutions
Swapping operator is applied in order to find the optimum solution. As this operator, the sequences of two parts, which are selected randomly, are being exchanged on each machine. For instance, suppose that the sequence of six parts on four machines is given as in Fig. 4. Exchanging the sequence of parts 5 and 6 on the first machine implies the modified solution as given by Fig. 5.

Replacement of the solutions
Two conditions are set to replace the new modified solution with the initial one: (a) if the modified solution improves the objective function value, it will be replaced with the initial solution; (b) a fraction of nonimproving solutions are accepted in the hope of escaping local optimum while searching global optimum. The probability of accepting nonimproving solutions is calculated as the Cauchy function instead of that of the Boltzmann, which can reduce the probability of trapping in local optimum [16]. The probability of replacing the initial solution P is given by Eq. (37) noting that ΔF is equal to the gap between the objective function's values of the modified and the initial solutions and T is the temperature of the existing iteration.

Initial temperature and stopping rule
Some statistical tests have been developed in order to select a proper value for the initial temperature and cooling rate expressed in Section 7.1. The temperature decreases using the cooling rate so that the new temperatures are found through multiplying it by the current temperature. The algorithm stops if the current temperature reaches to a predetermined value.   6 Computational results

Tuning SA parameters
The cooling rate and the initial temperature are the two parameters of the presented heuristic to be tuned. Experimentally, we have considered three levels for these parameters. The values selected for each parameter are given as in Table 4.
Ten numerical problems are designed as given by Table 5. The given parameters are stochastic with given probability distribution functions. Each problem is run for 100 times.
Factorial tests with a=0.05 are implemented. The quality of the solution and the CPU run time are the two criteria, which are considered for evaluating the performance of the proposed heuristics. In each condition, the null hypothesis is the equality of the two comparing levels for the mentioned  criteria. If the P value <a, the null hypothesis is not accepted. Results are illustrated by Table 6. Based on the quality of the solution and the CPU run time, level three (0.985) and level one (100) are more appropriate for the cooling rate and the initial temperature, respectively.

Measuring the robustness
Twenty problems with varying sizes have been designed in order to measure the robustness of the proposed heuristics.
According to the number of parts (NP), number of machines (NM), and number of cells (NC), the numerical problems are divided into three categories: small, medium, and large sizes. Two types of SAs are implemented, which differ on their initial solution generation procedure. The initial solution is generated randomly for the first type while for the second type, the SPT method is used. The comparison among the solution obtained from CPLEX and the two types of SAs is given in Table 5. Each problem was run for 100 times. Z mean , which represents the average of the objective function values, are reported in Table 7 along with the average of the run times. Z best is the best value found for the objective function, and SD represents the standard deviation of the obtained solutions (from 100 runs). CPLEX is able to find the optimum solution for small-and medium-size problems. In small-size problems, the SA-based heuristics are capable of finding the optimum solution. SA heuristics also show a promising performance for the medium-size problems. There is an almost 3 % gap between the optimum solution and the best solution of the hybrid SA, while this gap for the non-hybrid SA is almost 8 %. Furthermore, there is a huge difference between the run time of CPLEX and SA heuristics.
In term of the quality of solution, hybridizing SA with the SPT algorithm can improve the quality of solution. The best solutions of the hybrid SA with SPT are more than 4 % better than the normal SA on average. The SA with SPT also preserves its superiority in term of standard deviation of the solutions. In terms of the CPU runtime, SA has an advantage over SA with SPT. Figures 6, 7, 8, and 9 depict comparisons between SA and hybrid SA with SPT performances.

Conclusion and future research
Cell scheduling problem with routing flexibility and handling times was considered in this paper. A mathematical model was developed to minimize the makespan with sequencedependent family setup times. The optimum solutions are found even for bigger sizes of the problem converting the initial nonlinear model to a linear form; furthermore, the results in term of CPU runtime are improved. Unfortunately, linearization cannot solve the problem of high CPU runtime for very large problems; therefore, based on the previous researches, we developed a SA-based heuristic in order to solve the problem. The performance of the initial heuristic can be improved hybridizing it with SPT algorithm. The results for a number of small-, medium-, and large-size problems indicate that the heuristics have a promising performance.
Improvement of the structure of the model from both optimization and applicability aspects can be a good opportunity for the future researches. We believe there is a possibility to define the variables and constraints in a different way that can smooth the optimization process. Different assumptions can be added to the model that are compatible with the realworld phenomena. More efficient algorithms can be developed to find better solution in a more reasonable time.