Simulation optimization using tabu search

Proceedings of the 2000 Winter Simulation Conference J. A. Joines, R. R. Barton, K. Kang, and P. A. Fishwick, eds. SIMULATION OPTIMIZATION USING TABU SEARCH Berna Dengiz Cigdem Alabas Department of Industrial Engineering Gazi University Celal Bayar Bulv. Maltepe 06570 Ankara/TURKEY ABSTRACT and Chiu (1993) modeled a manufacturing cell with a fixed flow pattern considering dynamic effects of machine Investigation of the performance and operation of complex breakdowns and job changes. Ozdemirel et al. (1996) systems in manufacturing or other environments, analytical proposed a knowledge-based system, namely design of models of these systems become very complicated. experiments for simulation, that assist an experienced Because of the complex stochastic characteristic of the analyst. A general review on simulation optimization was systems, simulation is used as a tool to analyze them. The given in Tekin and Sabuncuoglu (1998). trust of such simulation analysis usually is to determine the Recently, many studies using modern heuristic optimum combination of factors that effect the considered techniques for simulation optimization have been system performance. The purpose of this study is to use a encountered. Bulgak and Sanders (1988) used modified tabu search algorithm in conjunction with a simulation simulated annealing to optimize buffer sizes in automatic model of a JIT system to find the optimum number of assembly systems. Grefensette (1991) considered strategy kanbans. acquisition with genetic algorithms. Stuckman (1988, 1990) has discussed the use of a particular Bayesian global 1 INTRODUCTION search algorithm for optimizing a design via simulation. Stuckman (1991) has compared three classes of global When investigating the performance and operation of search algorithms such as genetic algorithms and complex systems in manufacturing or other environments, Bayesian/sampling algorithms for design optimization. analytical models of these systems become very Haddock et al. (1992) used simulated annealing to optimize complicated. Because of the complex stochastic parameter levels considering the total profit of an characteristic of the systems, simulation is used as a tool to automatic production system. Hall and Bowden (1997) analyze them. However, the major drawback of simulation presented a comparative study of direct search methods for practical applications is that is time consuming. such as tabu search (TS), evolution strategies, and the To overcome the limitations of simulation, Nelder Mead Simplex algorithm for simulation metamodeling was first proposed by Blanning (1975). optimization. Lutz et al. (1995) have built a simulation Major benefits of metamodeling were summarized by model of a manufacturing process and used TS, a heuristic Madu (1990). Yu and Popplewell (1994) presented a procedure, to optimize buffer location and storage size in general review of metamodels in manufacturing. Barton this manufacturing system. Dengiz et al. (2000) used a (1992) reviewed different methods for choosing a regression metamodel to optimize batch sizes in a Printed functional form for the metamodel relationship such as Circuit Board assemble line. polynomials, Taguchi models and generalized linear Due to the manufacturing companies interest, models. Madu and Kuei (1992) employed group screening researchers started investigating the JIT philosophy and and Taguchi models in the design of experiment stage of a also much work has been done to find number of kanbans multi-echelon maintenance float simulation. Regression required in a JIT system. In general, a JIT system, if metamodels for steady-state systems are used, and they can implemented properly, will result in increased productivity, also represent dynamic behavior in response to several reduced work-in-process (WIP), and higher product quality unexpected real-time events in manufacturing. Lin and depending on the environmental factors of the JIT system. Cochran (1990) analyzed the dynamic performance of a In JIT systems both the level of WIP and order lead-time hypothetical multi-station, multi-server assembly line. Lin are important performance parameters. Inventory control in 805  Dengiz and Alabas a JIT system is controlled by the number of kanbans 2 KANBAN-CONTROLLED SYSTEM allocated. Kanban, which is a card in Japanese is used to direct materials to workstations and passes information as The example concerns the manufacture of two products, to what and how much to produce (Wang and Wang 1991). which are labeled Part7 and Part8. The production system Kimura and Terada (1981) describe the operation of includes two workcenters which are treated as black boxes kanban systems and examine the accompanying inventory (i.e. no specific flow patterns and machines are defined fluctuations in a JIT environment. Rees et al. (1987) within the cells). The kanban controlled system considered studied empirical approaches for setting kanban levels in this study is shown in Figure 1. The products are dynamically. Bitran and Chang (1987) presented a manufactured in two adjacent workcenters. The first mathematical formulation of the kanban determination workcenter uses raw materials Part1, Part2, Part3, Part4 problem. The formulation assumes planning periods of and produces two intermediate products, Part5 and Part6. known length and finds the minimum feasible number of The second workcenter gets the Part5 and Part6 products kanbans. Deleersnyder et al. (1989) places the kanban from workcenter 1, and produces the end products Part7 determination problem into the context of the overall pull and Part8. Part5, Part6, Part7 and Part8 are produced in the system implementation problem. Monden (1981) described manner given in Figure 1. the model in equation 1 for setting the number of kanbans for the Toyota Motor Company. 1 2 3 4 τ i Di (1 + α ) ki = (1) work center 1 ni 5 6 where ki is the number of kanbans for part type i, ni is the container size, τi is the sum of the lead time, waiting Part5 = Part1 + Part2 PK5 PK6 Part6 = Part3 + Part4 time and kanban collecting time, and Di is the average demand rate for part type i. While the demand rate is known on average, some variability does exist due to the WK5 WK6 order sequence at final assembly and drifts in demand. Because the safety factor, α, in equation 1 is to handle variability, the problem is the selection of α. Askin et al. work center 2 (1993) proposed an economic approach for selecting ki, and accordingly α. Their objective is to minimize the sum 7 8 Part7 = Part5 + Part6 of inventory holding and backorder cost. They formulated Part8 = Part5 a continuous time, steady-state Markov model to determine the optimal number of kanbans to use for each part type at PK7 PK8 each workcenter in a JIT system. The model selects the proper safety factor in each case. Figure 1: The Kanban-Controlled Manufacturing System Fukukawa and Hong (1993) proposed a mixed integer programming approach to examine many factors which For each type of product at each output buffer, specific play an important role in determining the number of production kanbans defined as PK5, PK6, PK7, PK8, and kanbans in a JIT production system. Their objective withdrawal kanbans WK5, WK6 are used to signal requests function was to minimize inventory holding, outage and for part transfers between the workcenters. The interarrival miscellaneous operating costs. Muckstadt and Tayur rate of the customer orders, order quantity and time delays (1995) presented a heuristic approach to determine the in the system are given in Table 1. number of kanbans. Aytug et al. (1996) determined the number of kanbans in a pull production system using a Table 1: Input Distributions for the Problem regression metamodel. Hurrion (1997) found the Interarrival rate of the customer N(3, 0.6) approximate optimum number of kanbans by using neural orders network metamodel. Order quantity N(1,3) The aim of this study is to find optimum number of Part5 N(1, 0.2) kanbans in a JIT system using TS. The performance of TS Time Delays Part6 N(1, 0.2) algorithm is compared with the performance of a random search algorithm (RS) applied on the same problem. The Part7 N(2, 0.4) example chosen to demonstrate this approach is a stochastic Part8 N(1, 0.2) discrete-event system described by Aytug et al. (1996). 806  Dengiz and Alabas Problem assumptions are: where: • An infinite supply of parts 1 - 4 is in the first f: total cost under given combination of workcenter kanbans, • Transfer time for production kanbans (PKs) is wc: waiting cost per minute and per order, negligible tnop: total number of order processed, • Withdrawal kanban (WK) transfer times are also kc: kanban cost per minute, zero (the succeeding workcenter is located very maxcyc: upper limit on the average order cycle close) time, • Container size is one for all the parts cyc7: order cycle time of Part7 (obtained from • The kanbans are released as soon as their the simulation model), containers are empty cyc8: order cycle time of Part8 (obtained from • Customer orders (demand) are external the simulation model), act: average cycle time for both part times = • Demands for Part7 40% of the time and for Part8 0.4*cyc7 + 0.6*cyc8 60% of the time makespan: completion time of all the orders. • Machine capacities are fixed and set at a level so that they will not generate any bottlenecks within Total cost f is a function of the number of kanbans the investigated range of the input variables only. Order cycles cyc7 and cyc8 are also functions of the number of kanbans only. Makespan could also be As shown in Table 1, the normal distribution is used to determined as a function of the number of kanbans in the represent all time delays in the system to create a relatively system. The makespan values for different combinations stable environment, which is the one of the main were very similar with only a small standard deviation. assumptions of a kanban implementation. Thus, the makespan value is computed as an average value The objective function of this study is to minimize the experimentally and considered as a constant for each possibility of backorders among workcenters and to keep the replication. customer order cycle time at a reasonably low level. Order All definitions, input values and parameters of the cycle time is defined as the difference between the time of kanban controlled system defined above are considered as the completion of the order and the time of the arrival of the the same of the example given in Aytug et al. (1996). order. Using order cycle time as a measure of performance Under these conditions, our simulation model was built and summarizes the effect of all the internal factors in a kanban validated according to the simulation model results given system. The order cycle times for each end product in Aytug et al. (1996). constitute the response variable values, which are obtained Examining kanbans, PK5, PK6, PK7, PK8, WK5 and by a simulation model for Part7 and Part8 separately. WK6, for this example means that at least 11,664 kanban combination points may need to be studied if the results of 3 THE PROBLEM FORMULATION the simulation model have a small standard deviation (meaning that no simulation replication is necessary). If the A simulation model for the defined example was coded in problem is expanded to a more complex form such as PASCAL to obtain the average cycle time of Part7 and adding one additional workcenter at the end of the line, that Part8 for each kanban combination. The problem is given is, an additional four kanbans (WK7, WK8, PK9, PK10) in in equation 2. the range of 1 ≤ WK7 ≤ 3, 1 ≤ WK8 ≤ 3, 1 ≤ PK9 ≤ 3, 1 ≤ Minimize PK10 ≤ 3, then the problem would require examining over 944784 different kanban combination points. f = wc * act * tnop + kc * totkan * makespan Due to the combinatorial nature of the number of Subject to kanbans optimization problem, the problem can become 1 ≤ PK5 ≤ 6 highly complex as the number of kanbans increase. It is not 1 ≤ PK6 ≤ 3 practical from a computational point of view to search the 1 ≤ PK7 ≤ 6 (2) complete set of all combination points. The approach 1 ≤ PK8 ≤ 6 considered in this study avoids having to evaluate all kanban combinations by employing a TS algorithm. As 1 ≤ WK5 ≤ 6 mentioned previously, the developed TS algorithm 1 ≤ WK6 ≤ 3 employed interacts the simulation model of the JIT system. 0.4cyc7 + 0.6cyc8 ≤ maxcyc totkan = PK5 + PK6 + PK7 + PK8 + WK5 + WK6 PK5, PK6, PK7, PK8, WK5, WK6 are integers 807  Dengiz and Alabas 4 TABU SEARCH Then the best neighbor is selected as the new current solution, if the neighbor is not obtained via a tabu move. Tabu Search (TS) is a meta-heuristic that guides a local The tabu list drives the search to different regions of heuristic search strategy to explore the solution space the search space. After any increasing or decreasing move beyond local optimality. The local procedure is a search operation, the activated kanbans are recorded on the tabu that uses an operation called a move to define the list. Two different tabu lists called tabu_increase_start and neighborhood of any given solution. The neighborhood of tabu_decrease_start were built to record the tabu the current solution is explored and the best solution is conditions associated with the moves a selected increase or selected as the new current solution. The best solution in decrease move at iteration i, tabu_decrease_start the neighborhood is selected, even if it is worse than the (increase_move) = i, or tabu_increase_start current solution. This strategy allows the search to escape (decrease_move) = i. The tabu tenure of increasing moves from local optima and explore a larger fraction of the (ttim) and tabu tenure of decreasing moves (ttdm) are the search space. Tabulist includes recently selected solutions number of iterations that forbid a kanban number to be that are forbidden to prevent cycling in the search process. increased and decreased, respectively. Any tested increase If the move is present in the tabulist, it is accepted only if it or decrease move is tabu at current iteration j (j > i), if decreases the objective function value below the minimal tabu_increase_start(test_increase_move) + ttim ≥ j or level so far achieved according to an aspiration level. It tabu_decrease_start(test_decrease_move) + ttdm ≥ j. This was originally proposed by Glover (1989, 1990). structure is also called short term memory, which is recency based (Glover 1989). An aspiration criterion was 4.1 Implementation of the TS Algorithm used to decide when the tabu rule can be overridden. The aspiration criterion used in this study removes the tabu The notations used in the developed algorithm are condition when any tested move yields a better solution introduced below: than best solution obtained so far. The parameters values of TS algorithm were y0: initial solution determined experimentally to be ttim = 3 iterations and y: current solution ttdm = 4 iterations. The developed TS algorithm is y′ : neighbor solution terminated when a chosen maximum iteration number is y′best: best neighbor solution reached. The steps of the TS algorithm are follows: ybest: best solution M(y): a move that yields solution y Algorithm: ttim: tabu tenure of increase moves Step 1. Choose the initial solution y0. Current ttdm: tabu tenure of decrease moves solution y = y0 and best solution ybest = y0. i = 0 and start with empty tabu lists. Solutions of the problem (candidate kanban Step 2. Repeat combinations) are represented by an array with six elements. This six elements include the number of kanban Step 2.1. Generate the neighbors, y′, for current solution y and call the for each kanban type, respectively, PK5, PK6, PK7, PK8, simulation model to calculate the WK5, WK6. Increasing and decreasing feasible moves were used to obtain neighbors of any solution. The possible cost function f(y′). neighbors of a sample solution [2, 3, 4, 3, 5, 3] are given in Step 2.2. Select the best neighbor y′best. If Table 2. f(y′best) < f(ybest) then ybest = y′best and go to step 2.4. Table 2: The Neighbors of Solution [2, 3, 4, 3, 5, 3] Step 2.3. If (M(y′best) is tabu) and (f(y′best) > [1, 3, 4, 3, 5, 3] [2, 3, 3, 3, 5, 3] [2, 3, 4, 3, 4, 3] f(ybest)) then f(y′best) = ∞ and go to step 2.2. else current solution y = [3, 3, 4, 3, 5, 3] [2, 3, 5, 3, 5, 3] [2, 3, 4, 3, 6, 3] y′best. [2, 2, 4, 3, 5, 3] [2, 3, 4, 2, 5, 3] [2, 3, 4, 3, 5, 2] Step 2.4. i = i + 1. Keep M(y′best) on [2, 4, 4, 3, 5, 3] [2, 3, 4, 4, 5, 3] [2, 3, 4, 3, 5, 4] associated tabu list for associated tabu tenure. The initial solution in the TS algorithm was randomly Step 3. Until i ≥ imax. selected and starting from the initial solution, all possible neighbors of the current solution are examined at each 5 RESULTS AND ANALYSIS iteration, because the number of neighbors is not too large. The TS algorithm calls the simulation model to compute In this study, the performance of TS algorithm is compared the total cost correspond to considered neighbor solution. with the performance of RS algorithm. RS is the simplest 808  Dengiz and Alabas heuristic search method that just samples solution space 6 CONCLUSIONS randomly. Thus, randomly generated kanban combination is used as input value for the simulation model of kanban - In this study, a simulation searched heuristic procedure controlled system to obtain objective function value. The based on TS was developed and compared with a RS algor- objection function value is compaired with the minumum ithm applied on the same problem considering both objective function value yet encountered. If it is smaller solution quality and computational efficiency for determin- than the the best one, it is stored as the new best solution. ing the optimum number of kanbans to meet production de- This process is repeated until a predetermined number of mands in a JIT system. Results indicate that TS algorithm solutions have been generated. outperforms the RS algortihm searching only 0.249% of The results show that the TS algorithm and the RS the solution space of this problem. The result encourages algorithm find the same kanban combination [1, 1, 3, 3, 1, 1] us to use TS method for simulation optimization. as the best result, while searching 29 and 794 solutions, respectively. The comparison among the heuristics is REFERENCES presented in Table 3 on the basis of three values: coefficient of variation of cost that is found by the relevant heuristic in each Askin, R.G., M.G. Mitwasi, and J.B. Goldberg. 1993. replication, number of solutions searched by the relevant Determining the number of kanbans in multi-item just- heuristic until finding the best solution (best kanban in-time system. IIE Transactions, 25(1): 89-98. combination that gives minimum cost) and fraction searched Aytug,H., C.A. Dogan, and G. Bezmez. 1996. Determining %. In order to ensure a comparable computational effort the number of kanbans: a simulation metamodeling devoted to each heuristic, the stopping criterion has been approach. Simulation, 67( 1): 23-32. defined as a number of solutions visited. This number has Barton, R.R. 1992. Metamodels for simulation input- been set to 900 in this study. Figure 2 represents the output relations. in Proceedings of the 1992 Winter convergence of the TS and the RS algorithms. This figure and Simulation Conference, ed., J.J. Swain, D. Goldsman, also Table 3 showed that the TS algorithm performed better R.C. Crain, and J.R. Wilson, 289-299. than the RS with the same number of solution visited, TS Bitran, G.R., and L. Chang. 1987. A mathematical converged very quickly visiting only 0.249% of the search programming approach to a deterministic kanban space. On the other hand, the nonlinear regression model used system. Management Science, 33(4): 427-441. by Aytug et al. (1996) to find optimum kanban numbers for Blanning, R.W. 1975. The construction and this problem fitted on the 64 simulation results (26 = 64) and implementation of metamodels. Simulation, 177-184. they also used 6 additional simulation runs to represent Bulgak, A.A., and J.L. Sanders. 1988. Integrated a validation of their regression metamodel. TS used in this modified simulated annealing algorithm with the study, only searched 29 points (means that 29 simulation runs) simulation of a manufacturing system to optimize to find optimum or very close optimum solution. buffer sizes in automatic assembly systems. In Proceedings of 1988 Winter Simulation Conference, Table 3: Results of the Computational Experiments ed., M. Abrams, P. Haigh, and J. Comfort, 684-690 TS RS Deleersnyder, J.L., T.J. Hodgson, H. Muller, and P.J. Coefficient of variation 0 11.8 O’Grady. 1989. Kanban controlled pull systems: an Solution searched 29 798 analytical approach. Management Science, 35(9): Fraction searched % 0.249 6.842 1079-1091. Dengiz, B., and K.S. Akbay. 2000. Computer simulation of a PCB production line: metamodeling approach. 450.000 International Journal of Production Economics, 63: 195-205. 425.000 Fukukawa, T., and S-C. Hong. 1993. The determination of 400.000 the optimum number of kanbans in a just-in-time production system. Computers and Industrial Cost 375.000 Engineering, 24(4): 551-559. 350.000 Glover, F. 1989. Tabu search Part I. ORSA Journal on Computing, : 190-206. 325.000 Glover, F. 1990. Tabu search Part II. ORSAJournal on Computing, 2: 4-32. 300.000 Grefensette, J.J. 1991. Strategy acquisition with genetic 0 60 120 180 240 300 360 420 480 540 600 660 720 780 840 900 Algorithms. In Handbook of Genetic Algorithms, ed., Solution Searched TS RS L. Davis, 186-201. Figure 2: Convergence of the TS and RS Algorithms 809  Dengiz and Alabas Haddock, J., and J. Mittenthal. 1992. Simulation Stuckman, B.E. 1991. Comparison of global search optimization using simulated annealing. Computers methods for design optimization using simulation. In and Industrial Engineering, 22(4): 387-395. Proceedings of the 1991 Winter Simulation Hall, J.D., and R.O. Bowden. 1997. Simulation optimiza- Conference, eds., B. Nelson, W.d., Kelton, and G.M. tion by direct search: a comparative study. In Institute Clark, 937-943. of Industrial Engineers 6th Industrial Engineering Tekin, E., and I. Sabuncuoglu. 1998. Simulation Research Conference Proceedings, 298-303. optimization: a comprehensive review on theory and Hurrion, R.D. 1997. An example of simulation applications. Technical Report, Bilkent University, optimization using a neural network metamodel: Ankara, Turkey. finding the optimum number of kanbans in a Wang, H., and H. Wang. 1991. Optimum number of manufacturing system. Journal of Operational kanbans between two adjacent workstations in a JIT Research Society, 48: 1105-1112. system. International Journal of Production Kimura, O., and H. Tereda. 1981. Design and analysis of Economics, 22(3): 179-188. pull system, a method of multi-stage production Yu, B., and K. Popplewell. 1994. Metamodels in control. International Journal of Production Research, manufacturing: a review. Industrial Journal of 19(3): 241-253. Production Research, 32(4): 787-796. Lin, L., and J.K. Cochran. 1990. Estimating simulation metamodel parameters for unexpected shop floor real AUTHOR BIOGRAPHIES time events. Computers and Industrial Engineering, 19(1-4): 662-666. DENGİZ BERNA is a Professor of Industrial Engineering Lin, L., and F. Chiu. 1993. Manufacturing cell operating and Vice Dean of The Faculty of Engineering at the Gazi Characteristics. European Journal of Operational University. Her field of study is the modeling and Research, 69: 424-437. optimization of complex systems by simulation and/or Lutz, C.M. 1995. Determination of buffer location and size metaheuristics techniques such as genetic algorithms, tabu in production lines using tabu search. European search, and simulated annealing. Her research has been Journal of Operational Research, 106: 301-316. funded The Scientific and Technical Research Council of Madu, C.N. 1990. Simulation in manufacturing: a Turkey (TUBİTAK) , National Planning Center of Turkey regression metamodel approach. Computers and (DPT) and National Science Foundation (NSF). Dr. Dengiz Industrial Engineering, 18(3): 381-389. is a member of SCS-Society for Computer Simulations, Madu, C.N., and C. Kuei. 1992. Group screening and IEEE –Reliability Chapter, Operations Research Society of Taguchi design in the optimization of multi-echolen Turkey, Informatics Society of Turkey, and Statistic maintenance float simulation metamodels. Computers Society of Turkey. Her e-mail address is and Operation Research, 19(3): 95-105. <[email protected]>. Monden, Y. 1981. Toyota production system. Norcross, GA: IE&MPress ALABAS CİGDEM is a Research Assistant of Industrial Muckstadt, J.A., and S.R. Tayur. 1995. A comparison of Engineering. She received as a B.Sc. (1995), M.Sc. (1999) alternative kanban control mechanism: I. Background in Industrial Engineering from Gazi University. She is and structural results. IIE Transactions, 27(2): 140- currently enrolled Ph.D. program in Industrial Engineering 150. department at Gazi University. Her research concerns Ozdemirel, N.E., G.Y. Yurttas, and G. Koksal. 1996. modeling of complex systems using stochastic Computer aided planing and design of manufacturing optimization techniques. She is a member of the simulation experiments. Simulation, 67(3): 171-191. Operations Research Society of Turkey. Her e-mail address Rees, l.P., P.R. Philipoom, B.W. Taylor, and P.Y. Huang. is <[email protected]>. 1987. Dynamically adjusting the number of kanbans in a just-in-time production system using estimated values of lead-time. IIE Transactions, 19(2): 199-207. Stuckman, B.E. 1988. Design of general systems by methods of global search: a computer-aided engineering approach. In Proceedings of the 1988 AMSE International Istanbul Conference on Modeling and Simulation, 1B: 81-90. Stuckman, B.E. 1990. Design optimization using simulation and stochastic global search: a computer- aided engineering approach. Advances in Modeling and Simulation, 18(4): 13-33. 810