A dynamic programming – enhanced simulated annealing algorithm for solving bi-objective cell formation problem with duplicate machines

Article history: Received June 10, 2014 Accepted October 22, 2014 Available online October 22 2014 Cell formation process is one of the first and the most important steps in designing cellular manufacturing systems. It consists of identifying part families according to the similarities in the design, shape, and presses of parts and dedicating machines to each part family based on the operations required by the parts. In this study, a hybrid method based on a combination of simulated annealing algorithm and dynamic programming was developed to solve a biobjective cell formation problem with duplicate machines. In the proposed hybrid method, each solution was represented as a permutation of parts, which is created by simulated annealing algorithm, and dynamic programming was used to partition this permutation into part families and determine the number of machines in each cell such that the total dissimilarity between the parts and the total machine investment cost are minimized. The performance of the algorithm was evaluated by performing numerical experiments in different sizes. Our computational experiments indicated that the results were very encouraging in terms of computational time and solution quality. Growing Science Ltd. All rights reserved. 5 © 201


Introduction
In today's competitive market, manufacturing industries must be able to produce products with low production cost and high quality to deliver the products to customers on time.In addition, they should be able to respond quickly to changes in product design and demand with low cost.Traditional manufacturing systems, such as job shops and flow shops cannot provide the efficiency and flexibility simultaneously to adapt to such requirements (Kioon et al., 2009).As a result, cellular manufacturing (CM), an application of group technology (GT), has emerged to satisfy such production requirements in manufacturing systems producing medium-volume/medium-variety products.GT is a manufacturing philosophy that identifies and explores the similarities of product design and manufacturing process to decompose a manufacturing system into several subsystems for facilitating the production floor control.In fact, CM is a hybrid manufacturing system that joints the advantages of flow shops and job shops with characteristics such as reduced cycle times compared to jobs shops and increased flexibility and greater job satisfaction as compared with flow shops (Mungwattana, 2000).If a cellular manufacturing system (CMS) is designed properly, it leads to shorter cycle times, reduced material handling cost, reduced work-in-process inventories, reduced setup times, reduced tool requirements, better control over manufacturing processes and improved quality of products (Wemmerlöv & Hyer, 1989;Jeon & Leep, 2006).To implement an efficient CMS, several issues such as cell formation (CF), facility layout, production planning, scheduling, etc. should be considered in the design process.Among these issues, CF process is one of the first and most important steps in designing CMSs.It consists of identifying part families according to the similarities in the design, shape, and presses of parts and dedicating machines to each part family based on the operations required by the parts.This process is also addressed as the machine-part grouping problem (MPGP) in the literature.In the context of the CF problem, numerous approaches have been developed in the literature.Among those approaches, classification and coding systems, similarity coefficient-based methods, array-based clustering methods, mathematical programming, graph partitioning, artificial intelligence (AI) based methods and heuristic and meth-heuristic based algorithms are the frequently addressed approaches (Singh, 1993;Mungwattana, 2000).
As CF problems belong to the class of NP-hard combinational problems, most researchers have focused on implementation of meta-heuristic algorithms such as genetic algorithm (GA), simulated Annealing (SA), tabu search (TS), ant colony optimization (ACO), particle swarm optimization (PSO), etc.For instance, Boctor (1991) presented a 0-1 linear mathematical model for MPGP and aimed to minimize the number of exceptional elements (EEs).A SA algorithm was developed to solve the model.The proposed SA seems to able to find the optimal solution for 58 (64.4%) of the 90 solved problems.Chen and Srivastava (1994) formulated the CF problem as a quadratic programming model to maximize the total similarity between machines, subject to cell size limitation.They employed SA to solve the problem and concluded that the performance of the SA is better than the graph-partitioning heuristic.Balakrishnan and Jog (1995) proposed a parallel genetic TSP algorithm to minimize the number of EEs in the MPGP.In this study, the MPGP was represented as two TSPs by converting the similarity coefficients between the parts and between the machines into distances.Caux et al. (2000) addressed the problem of manufacturing cell formation with alternative processing routings, machine capacities, operation times and part demands.A hybrid approach combining the SA algorithm for the CF and a branch-and-bound (B&B) method for the routing selection was presented to minimize the number of intercellular movements.As the proposed algorithm uses the B&B method, it may not be efficient in solving large-sized problems or unconstrained problems both in terms of accuracy and computational time.Adil and Rajamani (2000) presented anon-linear mathematical model to investigate the trade-off between cell compactness and cell independence in terms of inter-cell and intra-cell move costs.A SA algorithm was used to solve the problem.Mungwattana (2000) studied the CF problem with alternative processing routings in dynamic and stochastic production requirements and utilized SA algorithm to solve the problem.Onwubolu and Mutingi (2001) developed a GA to solve the MPGP with two objectives, minimizing the total number of EEs and minimizing the total cell load-variation.The proposed GA appears to perform better than TSP-based heuristics.Solimanpur et al. (2004) presented a mathematical model for designing independent manufacturing cells with multiple objectives of minimizing the total dissimilarity, total processing cost, total processing time and total investment in the acquisition of machines.A GA with multiple fitness functions was used to find multiple solutions along the Pareto optimal frontier.Chiang and Lee (2004) developed a SA algorithm augmented with dynamic programming to solve the CF problem with an objective function of minimizing the inter-cell handling cost.Tavakkoli-Moghaddam et al. (2005) employed GA, SA and TS to solve a modified version of proposed problem by Mungwattana (2000).Their computational experiments indicated the superiority of SA over GA and TS.In addition, Tavakkoli-Moghaddam et al. (2008) modified the same work by considering reconfiguration and employed the SA algorithm to simultaneously minimize the inter-cell movements and machine costs.Their computational experiments showed that the gap between optimal and SA solutions is less than 4%.Arkat et al. (2007) developed a sequential CF model and solved it by SA and GA.They reported similar results for both methods.However, SA needed less computational time.(Saghafian & Akbari Jokar, 2009) extended the work of Chiang and Lee (2004) to intra-cell layout and developed a hybrid method based on SA, ACO and dynamic programming to minimize the total inter-cell and intra-cell handling cost.Banerjee and Das (2012) defined a new grouping efficacy in terms of a linear combination between the operation densities within the cells and the number of inter-cellular movements.A predator-prey GA was employed to solve the MPGP by considering this measure.Ghezavati and Saidi-Mehrabad (2011) applied queuing theory to formulate the CF problem with stochastic parameters.They assumed that each machine works as a server and each part is a customer and aimed to maximize the average probability that machines are busy in cells.A hybrid method based on combination of SA and GA was proposed to solve the problem efficiently.This hybrid method was compared against global solution obtained from B&B algorithm and a benchmark heuristic algorithm.They reported successful performance in any size of problems.Kao and Lin (2012) proposed a discrete PSO algorithm to minimize the number of EEs in the CF problem.They concluded that the PSO algorithm results in better solutions in comparison with the SA and TS-based algorithms.
In spite of many researches on the CF problem, most of the existing methods or algorithms in this area are not general enough to determine the number of machines in the cell design process.These approaches usually start with an initial machine-part matrix and then rearrange the columns and rows of this matrix to form part families and machine cells.In other words, those methods do not consider the capacity of machines and the processing information of parts such as demand and processing times in the cell design process.In this study, a bi-objective mathematical model in CMSs is presented to design independent manufacturing cells considering duplicate machines.The proposed model is to determine the part families and the number of machines in each cell such that the total dissimilarity between the parts and the total machine investment cost are minimized.Due to the complexity and NP-hard nature of CF problems, a hybrid solution approach based on combination of SA algorithm and dynamic programming is developed to solve the problem efficiently.After setting the parameters of the proposed algorithm, its performance is evaluated by performing numerical experiments in different sizes.

Model description
In our problem, parts are grouped into the part families based on the dissimilarity measure defined between them and considering the maximum number of parts allowed in a cell.In addition, machines are allocated to each part family such that the demand of all parts are satisfied and all the operations required by a part family are processed within a cell (i.e., the EEs are eliminated by duplicating machines).Two conflicting objective functions are considered in the problem.The first objective function, which minimizes the total dissimilarity between the parts are used to control the flow lineness of the system.From the other side, the second objective function is used to control the job shopness of the system by minimizing the total machine investment cost.In other words, the designer will be able to adjust the job shop-ness or flow shop-ness of the manufacturing system in terms of the weight of each objective in the problem.

Notations
The following notations are used in the formulation of the problem.

Sets:
, parts index , 1, … , ( is the number of parts) machines index 1, … , ( is the number of machine types) cells index 1, … , ( is the maximum number of cells allowed) Parameters: demand of part available time of machine type purchase price of machine type =1 if part visits machine ; 0 otherwise processing time of part on machine type dissimilarity coefficient between parts and maximum number of parts permissible in a cell (cell size limit) total dissimilarity between parts total machine investment cost , weights of and , respectively

Decision variables:
=1 if part is assigned to cell ; 0 otherwise number of machines of type assigned to cell

Problem formulation
In context of similarity and dissimilarity coefficients, extensive reviews can be found in (Yin & Yasuda, 2005;Yin & Yasuda, 2006;Garbie et al., 2008).Among various similarity coefficients, Jaccard similarity coefficient proposed by McAuley (1972) is the most stable and most often used similarity coefficient in the literature Yin and Yasuda (2005).It is also a very simple and effective measure in designing CMSs.In this study, the dissimilarity coefficient between two parts and is defined by subtracting Jaccard similarity coefficient from its upper bound of 1, and it is defined as follows: Now, the proposed bi-objective problem can be formulated as the following non-linear integer model: min , (2.1) min , .
(2.2) Subject to: Objective function (1.1) minimizes the total dissimilarity and objective function (2.2) minimizes the total machine investment cost.Constraint (3) ensures that each part is assigned to one cell.Constraint (4) ensures that the cell size limit is not violated.Constraint (5) determines the number of each machine type in each cell.Lastly, constraints ( 6) and ( 7) indicate the type of decision variables.
The mathematical model presented above is non-linear, due to the presence of the product term in objective function (2.1).Most non-linear problems are usually much harder to solve optimally than linear problems.Here a method from Kaufmann and Broeckx (1978) which has the smallest number of variables and constraints is applied to linearize the model.To do this, objective function (2.1) is rearranged as follows: ∑ ∑ , .Now, by introducing a new set of variables to replace the ∑ product terms, and also defining constraints ( 9) and ( 10), objective function (2.1) is linearized as follows: min , .
(8) Subject to: To convert the proposed bi-objective model into a single objective one, a weighted sum of objective is used as follows: Subject to: (3)-( 7), ( 9) and ( 10).
In objective function (11), parameters ′ and ′ are defined as follows: where L and U respectively denote the lower and upper bounds of , and L and U respectively shows the lower and upper bounds of .To determine L and U we can solve the problem with the following objective function: min (Where is a small enough number).Also, U and I can be obtained by solving the problem with the following objective function: min .

Proposed enhanced SA algorithm
SA is a stochastic search method, which imitates the physical annealing of solid for solving combinatorial optimization problems.It has been known that the CF problem is one of the NP-hard combinational problems (Kazerooni et al., 1997).It means that obtaining optimal solution in a reasonable computational time is difficult, especially for large-scale problems.In recent years, SA algorithms have been successfully applied by researchers for solving the CF problems (Chiang & Lee, 2004;Tavakkoli-Moghaddam et al., 2005;Arkat et al., 2007;Ghosh et al., 2011).These have motivated us to employ the SA algorithm for solving the proposed problem.In addition, a dynamic programming algorithm is developed to improve the performance of the SA.The main elements of the proposed solution approach are explained below.

Generating initial solution
In the SA implementation, each solution must be represented by a coding scheme.In this work, a string of randomly generated integer values is used to represent a solution.If a problem involves parts, a permutation of integer values with the length of bits is needed to encode the solution.This string is associated with the sequence of parts at the machine-part matrix, which must be partitioned into part families and machine cells.In this way, a dynamic programming algorithm is developed to solve the partitioning problem.It should be noted that applying this coding scheme yields only feasible solutions and also leads to the reduction of the string length needed to represent a solution Chiang and Lee (2004).The proposed solution approach has been illustrated in Fig. 1.

Evaluating solutions by dynamic programming
The proposed scheme for encoding the solutions gives only the permutation of parts in the machinepart matrix and not the part families and the number of machines required in each cells.As it has been assumed that the machine cells must be formed independently (i.e., no material movement is allowed between the cells), the partition problem can be solved by dynamic programming algorithm.
Here, a dynamic programming algorithm is presented to partition the permutation of parts into the part families and determine the number of machines in each cell under the cell size limit ( ) and the maximum number of cells ( ) constraints.

Dynamic programming algorithm
In this section, we show how a dynamic programming approach can be used to solve the partition problem.Let denotes a permutation of parts.Now, the partition problem can be stated as follows.
In a permutation of parts, a breaking node is the part, which forms a part family (or one cell) and the node after the breaking node is the beginning for forming the next part family.Let be the order index of node in permutation , and denotes the part index which has placed in order .The partition problem is to find a set of breaking nodes ( ) which partition the permutation of parts into part families in such a way that the weighted sum of the total dissimilarity and the total (13) subject to: where , calculates the increased dissimilarity from breaking node 1 to breaking node , and , calculates the increased machine investment cost from breaking node 1 to breaking node .These parameters are calculated as follows: In Eq. ( 19) the symbol indicates the smallest integer value bigger than .
Objective function ( 13) minimizes the weighted sum of the total dissimilarity between the parts and the total machine investment cost.Constraint ( 14) ensures that each part family must contain at least one part.Constraint (15) makes sure that the cell size limit is not violated.Constraint (16) restricts the number of cells allowed to be formed.Finally, constraint (17) represents that the first cell must be partitioned from the first breaking node.
By using the dynamic programming, the partition problem can be sequentially solved in stages from 1 to .Let be the index of stages, , indicates the objective function value at stage , when at this stage the permutation of parts is partitioned from node 1 to , and * be the optimum objective function value at stage 1, when its breaking node is .Therefore, the partition problem at stage can be stated as follows: where 0, * 0 and * min , .In addition, constraints ( 21) and ( 22) control the cell size limit and the maximum number of cells constraints to avoid infeasible solution.
Note that, the proposed dynamic programming approach partitions a string of solution into exactly cells.Therefore, the optimum objective function value, * , and the optimum number of cells, * , for permutation are obtained as follows: The pseudo code of the proposed dynamic programming algorithm is given in Fig. 2.
Fig. 2. Pseudo code of the proposed dynamic programming algorithm for evaluating the optimal objective function value of a permutation of parts

An illustration
To illustrate the implementation of the dynamic programming algorithm, an example with 10 parts and 10 machines is solved.The input parameters for this example are given in Table 1.For simplicity, the purchase prices of all machines are assumed to be 1$.Also, the design factors are assumed as follows: 4, 4, ′ 0.2 and ′ 0.8.The proposed dynamic programming algorithm is implemented on typical solution 10, 5, 1, 7, 9, 3, 2, 6, 4, 8 to partition it into part families and determine the number of machines in each cell.For this permutation, the increased dissimilarity, , , and the increased machine investment cost, , of partitioning a cell from breaking node to breaking node were calculated and reported in Table 2.As the maximum number of cells is equal to 4 cells (i.e., 4), the dynamic programming is solved in 4 stages.The summary of the computations are reported in Tables 3-6.According to the results of Tables 5 and 6 we conclude that partitioning permutation 10, 5, 1, 7, 9, 3, 2, 6, 4, 8 into 3 cells results in better objective function value, because * 10 * 10 .Therefore, the optimal objective function value, * , and the optimal number of cells, * , are obtained as 15.33 and 3, respectively.The final solution is shown in Table 7.

Cooling schedule and moving to a neighboring solution
SA algorithm works with a controlled cooling schedule, which is also called the annealing schedule.
Starting from initial temperature , the temperature is gradually decreased through a pre-determined cooling schedule.The most commonly used cooling function is the geometric decrement function, , first proposed by Kirkpatrick et al. (1983).Where is the cooling rate and takes value between 0 and 1.In this study, this function is used to decrease the temperature.The value of the temperature at the beginning of the schedule should be large enough so that most of the initial movements are accepted.However, as the temperature reduces, the probability of accepting a nonimproving solution reduces.To calculate the initial temperature, we modify and combine the approaches proposed by (Tavakkoli-Moghaddam et al., 2008;Ghezavati & Saidi-Mehrabad, 2011).So, the initial temperature can be calculated as follows: where and are two randomly generated solutions at th trial.
At each iteration (temperature), a generation mechanism called Move is applied to transform the current solution into a neighboring (new) solution.Three move operators are proposed, namely Swap, Change and Inverse operators.The Swap operator swaps the order of two randomly selected parts.
The Change operator changes the order of a randomly selected part.The Inverse operator reverses the order of parts between two randomly selected points.To illustrate these operators assume that the current solution is Current 5, 3, 6, 1, 8, 2, 4, 7 , if the order of part 2 is changed to 3 (Change operator), the neighboring solution becomes New 5, 3, 2, 6, 1, 8, 4, 7 .If parts 3 and 8 are chosen for swapping (Swap operator), the neighboring solution becomes New 5, 8, 6, 1, 3, 2, 4, 7 .Lastly, if the orders between parts 1 and 4 are inverted (Inverse operator), the neighboring solution becomes New 5, 3, 6, 4, 2, 8, 1, 7 .It should be noted that these operators are performed on the solution independently to derive a neighboring solution.

Set
1: New ← Change Current ---2: New ← Inverse Current ---3: New ← Swap and Change Current ---4: New ← Swap and Inverse Current ---5: New ← Change and Inverse Current ---6: New ← Swap, Change and Inverse Current --End case Current .If the change in each transition represents a reduction in the objective function value i.e., Δ 0, the transition to the new solution is accepted.Otherwise, the nonimproving solution is accepted with a specified probability function exp Δ ⁄ .Accepting nonimproving solutions enables the SA to avoid the entrapment in local optima.This mechanism at each temperature is repeated until accepted transitions are met.Where is proportional to the length of string used in the solution encoding (i.e., ).

Stopping criteria
The proposed SA algorithm is terminated when a specified number of iterations, , is reached.The detailed steps in the implementation of the proposed SA algorithm are given in Fig. 3.

Experimental studies
The proposed dynamic programming-enhanced SA algorithm was coded in the Embarcadero Delphi XE programming language and implemented on a Windows 7 PC with 2.4 GHz CPU and 2 GB RAM.As the quality of solution is usually sensitive to the SA parameters, experimental tests are conducted to set the SA parameters.Then, the performance of the SA is evaluated in two sections.In the first section, a set of randomly generated instances are solved by the SA, and the results are compared with the solutions obtained by the B&B method.In the second section, numerical examples adopted from the literature are solved by the SA algorithm and the results are compared with those from the literature.

Parameter settings for SA
As widely known, settings of SA parameters critically affect the solution efficiency and effectiveness.
The parameters of the proposed SA algorithm include: cooling rate ( ), number of accepted transitions at each iteration ( ) and number of iterations ( ).The value of each parameter is selected from a set of predefined values given in Table 8.It should be noted that these values have been obtained based on our initial computational experience and the result of similar algorithms in the literature.Five test problems with 15, 20, 25, 30, 35 are used in the parameter settings.The input parameters of these test problems are generated randomly according to Table 9.In this table, the term " " implies uniform distribution and the term '∑ ' denotes the number of operations of each part which is selected randomly between 3 and 5.All the problems are investigated for 0.5.To calculate ′ and ′ , we first determine L , U , L and U .In this way, the SA algorithm is solved two times regarding the following objective functions: min and min (the first objective function gives U and L , and the second one gives U and L ).We then obtain ′ and ′ using Eq. ( 12).For each combination of parameters, each example was solved 10 times (in total 5 10 5 3 2250 experiments were done) and Minitab 16 software was used to analyze the results.The effects of parameters on the solution quality and computation time are plotted in Figs. 3 and 4, respectively.To select the best combination of parameters, we consider the trade-off between the solution quality and computation time.As it can be seen in Figs. 4 and 5, decreasing from 0.7 to 0.6, increasing from 30 to 40 and increasing from 3 to 4 does not improve the solution quality meaningfully, but it significantly increases the computation time.According to these plots, the parameters of the SA algorithm are obtained as follows: 0.7, 30 and 3 (the trade-off points have been highlighted in Figs. 4 and 5).9.For all the problems we assume that 0.5, and calculate ′ and ′ as explained in Section 4.1.As it was mentioned before, the CF problems are NP-hard.Hence, some problems may not be solved optimally by the Lingo in an acceptable amount of time.Thus, the solver is interrupted after 7200 seconds (2 hours) and the optimality gap (the relative gap between current feasible solution and best bound on optimal solution) is reported.From the other side, due to the stochastic nature of the SA algorithm, each problem is solved 30 times and the best result is considered for comparison.Table 10 demonstrates the results of comparison.The columns of this table are defined as follows: the best objective function value obtained by the Lingo ( B&B ), the bound on the objective function value ( Bound ), the time spent on solving the problems by the Lingo ( B&B ), the relative gap between the objective function value and its bound (Opt.Gap), where Opt.Gap The problems with 10, 11, … , 15 were solved optimally by the Lingo.However, the remaining problems were not solved optimally in 2 hours.As it can be seen in column "Opt.Gap", by increasing the problem size the optimality gap increases.This implies the complexity of the problem.From the other side, the results show that in all the problems the solution of the SA is better than or at least the same as that obtained by the Lingo.For these problems, the average times that the SA algorithm found better solution than B&B is almost equal to28 out of 30 runs.In addition, the standard deviation of the objective function values are very small and the mean of the objectives function values are very close to the best objective function values.These demonstrate that the SA algorithm is able to consistently produce good solutions.Moreover, the results indicate that the SA algorithm is able to solve the problems in a very short time (less than 4 seconds) compared to the B&B algorithm.

Comparison with other studies
In this section, twelve problems in different sizes have been selected from the literature.For some problems all of the necessary information was not available in the source papers, so we added the required parameters to the original problems.Table 11 shows the source, size and the parameters added to each problem.For all the problems we assume that 0.5, and calculate ′ and ′ as explained in Section 4.1.These problems are solved by the proposed SA algorithm and the results are compared with those reported in the literature.As it was mentioned in Section 2, most of studies dealing with the CF problem do not consider machine duplication in the problem formulation, hence to be able to compare the results we have to calculate the number of each machine type in each cell regarding the part family results available in the literature.The summary of comparison is given in Table 12.In this table column "Imp.(%)" shows the improvement percent in the objective function value and is defined as follows: Imp  The results of Table 12 show that in all the problems (except for problems 3, 5 and 9) both the total dissimilarity, SA , and the total machine investment cost, SA , obtained by the SA algorithm are better than those obtained for the solution reported in the literature.In addition, it can be seen that in all the problems the objective function value of the SA algorithm, SA , is considerably better than that calculated for the solution reported in the literature.For these problems, the average improvement in the objective function value is equal to 13.47%.

Conclusions
In this paper, a SA algorithm was developed to solve a bi-objective CF problem with duplicate machines.In the proposed SA algorithm each solution was represented by a permutation of parts and a dynamic programming algorithm was employed to partition this permutation into part families and determine the number of each machine type in each cell such that the total dissimilarity between the parts and the total machine investment cost are minimized.After setting the SA parameters using statistical experiments, a set of randomly generated instances were used to compare the SA algorithm against the B&B algorithm.The results indicated that the SA algorithm is able to find optimal solution in a very short computational time.In addition, several numerical examples adopted from the literature were solved by the proposed SA algorithm and the results were compared to those from the literature.It was showed that the proposed SA algorithm gives much better results than the others for all problems with an average improvement of 13.47% in the objective function value.
Partitioning the permutation of parts into the part families and determining the number of each machine type in each cell using dynamic programming Design factors: 19 10 13 19 16 16 12 machine investment cost are minimized.Based on these definitions, the partition problem can be formulated as the following integer programming model:

Fig. 3 .
Fig. 3. Pseudo code of the proposed SA algorithm

Fig. 4 .
Fig. 4. Effect of parameters on solution quality

Table 1
Data of the example used in the illustration of the dynamic programming algorithm

Table 8
Experimental factors and their levels

Table 9
Data set generation Best ), the mean of CPU times in 30 runs of the SA ( SA ), the mean of standard deviation of the objective function values ( ), the number of times that the solution of the SA is better than that obtained by the Lingo ( ), and the relative gap between Best and B&B (Gap), where Gap where SA is the best objective function value found in 30 runs of the SA algorithm and Lit is the objective function value calculated for the solution reported in the literature.

Table 11
Characteristic of the problems selected from the literature

Table 12
Result of Comparison for the problems selected from the literature