A Metaheuristic Based Self-Reasoning System for Assembly Sequence Automation in CIM

©BEYKENT UNIVERSITY

A METAHEURISTIC BASED SELF-REASONING

SYSTEM FOR ASSEMBLY SEQUENCE

AUTOMATION IN CIM

Cevdet GOLOGLU

Division of Design and Machine Building, Faculty of Technical Education Karabuk University, 78050 Karabuk

E-mail: cgologlu@hotmail.com

Received: 6 February 2007, Revised: 3 July 2007, Accepted: 27 November 2007

ABSTRACT

Assembly sequencing task is a time consuming process. In traditional manufacturing, it requires human planning in order to decide which sequence will be traced to assembly the parts. Even in some modern Computer Integrated Manufacturing (CIM) environments coupled with assembly robots after human operated sequence determination, the route of assembly robot are coded to robot control unit. Automation of this process, of course, will eliminate the error caused by human factor and drop the cost of product. Thus, the objective of paper is to determine intelligently configuration of robot movements for a number of parts that will be assembled. A prototype of gantry-type robot operating as a material handing device for an assembly automation system was designed and produced. Three metaheuristic techniques, Genetic Algorithm (GA), Simulated Annealing (SA) and Tabu Search (TS), were employed to decide the shortest path of assembly sequence. Then, the optimum assembly configuration found by one or more techniques was passed automatically to a control unit to carry out assembly tasks. The algorithms are shown to be effective regarding the results obtained. The model developed can be employed as a self-reasoning system for facilitating automation of assembly in a CIM environment.

Key Words: Genetic Algorithm (GA), Simulated Annealing (SA), Tabu Search

BİLGİSAYAR TÜMLEŞİK İMALATTA MONTAJ

OTOMASYON SIRALAMASI İÇİN BİR

METAHURİSTİK TABANLI OTOMUHAKEME

SİSTEMİ

ÖZET

Montaj sıralama görevi zaman alıcı bir süreçtir. Geleneksel imalatta parçaların montajında hangi sıranın izleneceğinin kararlaştırılması için insan planlamasına gereksinim vardır. Hatta bazı montaj robotları ile donatılmış modern Bilgisayar Tümleşikli İmalat (BTÎ) ortamlarında, montaj robotunun rotası insan tarafından montaj sırasının belirlenmesinden sonra robotun denetim birimine kodlanır. Bu süreçteki otomasyon, insan etmeni nedeniyle oluşan hataları kesinlikle azaltacak ve ürün maliyetini düşürecektir.

Bu makaledeki hedef, montaj edilecek parçalar için robot hareketlerinin konfigürasyonunu zeki bir şekilde belirlemektir. Bu amaçla otomatik montaj sistemi için malzeme taşıma-konumlama aygıtı olarak kullanılacak vinç tipi bir prototip robot tasarlanmış ve üretilmiştir. Üç adet metaheuristik teknik: Genetik Algoritma (GA), Tavlama Benzeşimi (TB) ve Yasak Araştırması (YA), en kısa montaj sıralama yoluna karar vermek için kullanılmıştır. Sonra, bir yada daha fazla metaheuristik teknik ile bulunan optimum montaj konfigürasyonu otomatik olarak montaj işlerini yerine getirmek için robot denetim birimine gönderilmiştir. Elde edilen sonuçlar ışığında algoritmaların etkili olduğu göstermiştir. Geliştirilen model bir BTİ ortamındaki montaj otomasyonunu kolaylaştırmak için kendi kendine muhakeme sistemi olarak kullanılabilir.

Anahtar Kelimeler: Genetik Algoritma (GA), Tavlama Benzeşimi (TB), Yasak

Araştırması (YA), Vinç tipi robot, Montaj otomasyonu.

1. INTRODUCTION

Production machines add value to the products by processing raw materials. Many manufactured parts are required to fit into a complete machine or structure. The labor cost for assembly varies between 50% and 70% of the total labor cost for manufacturing the product [1]. One of the purposes of CIM is to automate the associated processes. Many of different type robot configurations are employed in the processes. Tooling attached to Z-axes, the robots with cartesian configuration offer the advantage of rigidity to assembly operations. The enclosed rectangular work envelope of this type is often used

to move parts from conveyor systems into production machines and also used to assembly parts in assembly stations.

Assembly sequencing task is a time consuming process. In traditional manufacturing, it requires human planning in order to decide which sequence will be traced to assembly the parts. Even in some modern environments coupled with assembly robots after human operated sequence determination, the route of assembly robot are coded to robot control unit [2]. Automation of this process, of course, will eliminate the error caused by human factor and drop the cost of product. Researchers [3-6] have studied how sequences of assembly of parts can be presented. Emerging software engineering tools have been applied successfully to solving optimization problems in the presence of many local extrema. The problem of assembly sequence generation, which is kind of TSP, requires finding acceptable approximations like of many other NP complete problems.

This paper is structured as follows. The next section gives brief explanation of the gantry-type robot produced. In Section 3, how the system intelligently carries out the task, the algorithms used, and the problem formulation are explained in details. The performance study and evaluation of GA, SA and TS on a sample short path finding of an assembly automation process are discussed in Section 4 while the last section concludes the paper.

2. DESIGN OF A GANTRY-TYPE ROBOT

Many manufacturing and assembly operations require the partially finished parts or raw materials being oriented, inserted and moved from one location to another. As a result, some auxiliary production hardware such as material handling systems, robot arms and storage systems are needed. Among them robot manipulators have significant importance in manufacturing automation. The basic mechanical configurations of robot manipulators can be categorized as cartesian, cylindrical, spherical, and articulated [7].

A gantry-type robot within cartesian geometry, which has a rectangle work envelope moving in the X, Y and Z directions, was designed and produced [8]. It can move its gripper to any position within its working volume. The actual working volume of the robot is less than the theoretical volume of the robot resulting from mechanical design constraints. Figure 1 shows the prototype gantry-type robot produced for proposed implementation.

Figure 1. Prototype gantry-type robot produced and system environment

3. DETERMINING OF ASSEMBLE SEQUENCE

In assembly sequencing, a number of parts should be assembled such a sequence that a shortest assembly route can be achieved. The size of the problem affects the selection of optimization technique. When it is sufficiently large, sub-optimal techniques, stochastic or heuristic rules, that avoid exploring the whole state space, must be introduced to search for solutions [9, 10]. Metaheuristic algorithms are able to explore efficiently a large number of possible solutions [11]. A self-reasoning system, employing GA, SA and TS algorithms, was integrated with the prototype to find the shortest path of assembly sequence. The system acquires the input data, coordinates of the parts to be assembled, from computer database (Figure 2). Then, the metaheuristic algorithms, GA, SA and TS, all run to find the shortest path of assembly sequence for each time a problem is encountered.

The system, finally, takes the optimum (the shortest path resulted by the algorithms) assembly configuration found one or more algorithms to pass automatically to the control unit to carry out assembly tasks. Since all the methods are stochastic and they may not guarantee of optimality, however, its results have been compared with the existing solutions known.

Metaheuristic Techniques

Genetic Algorithm (GA) Simulated Annealing (SA)

Tabu Search (TS)

•

Figure 2. Schematic view of the self-reasoning system.

3.1. Genetic algorithm

GA, originally developed by Holland [12], is a robust parallel search technique and models natural selection and the process of evolution. The algorithm uses natural selection mechanism in evolving individuals, which are chromosomes usually represented as bit strings, through selection and recombination operators. Evolution is conducted from one population of chromosomes (or individuals) to a new population in keeping the population size constant throughout a number of iteration. Evaluation of each string is based on a problem dependent fitness function that helps determine the survivability in natural selection. Fitter individuals have higher probability of being selected. The chromosomes selected form the new population. Then, a number of genetic operators are applied to the chromosomes to generate offsprings. Figure 3 depicts outline of the genetic algorithm.

In recombination, the genetic operators, which are crossover and mutation, generate new individuals by combining the features of selected parents. Crossover, applied with high probability, allows information exchange between the selected pair of chromosomes by choosing random points and exchanging the substrings set by the points. Mutation, applied with very low probability, flips a randomly chosen bit. This survival of the fittest process continues until a termination criterion is satisfied. The best fit states optimal or near optimal result.

1. Generate random initial population P of solutions (n chromosomes) 2. Evaluate the fitness f(x) of each chromosome x in P, find the Best 3. REPEAT

3.1 Select two parent chromosomes

3.2 Cross over the parents with a crossover probability pc

3.3 Mutate new offspring with a mutation probability pm

3.4 Place new offspring in a new population UNTIL new population is complete

4. Replace new generated population for a further run 5. IF Stopping criterion = False THEN Go to 2

6. ELSE Return to Best solution.

Figure 3. Outline of the genetic algorithm.

3.1.1. Reproduction

GA implementation of the problem is as follows. GA usually uses binary coding for representation of a possible candidate solution. However, if it is used in ordering and for complex problems, a chromosome requires a longer gene chains resulting big computation time in decoding the chromosome. Thus, the use of decimal code is more appropriate for solving sequencing problems. In the study decimal code has been employed to build chromosomes. In other words, the possible solutions are permutations and represented as strings of integers. The set of chromosomes which comprise the initial population is generated at random. The population size in G A should be considered carefully. If the population is too small, it is very likely that premature convergence of the population to a local optimum occurs. Thus, the population size of a couple of hundreds in the problem is enough to escape from premature convergence. The fitness function in the problem is the total minimum distance between the assembly points to be visited.

Fmm = ( 1 )

i=1

where n is number of assembly points with the convention that n+1 is identified with point 1 and each point has coordinates of (xh yi).

Selection for reproduction in GA is done in several ways. The simplest scheme is construction of a roulette wheel with slots sized according to fitness values. A random number is generated and the chromosome whose interval sweeps the number is selected. The spinning process continues until all members of the new population are assigned. The new population carries individuals from old

population dominating the stronger strings and giving fewer chances to the weaker ones. However, there are some disadvantages as the selection process can lead to premature convergence on local optimum. It is due to the outstanding individuals will, in general, be reproduced more frequently while some of the poorer strings will probably not be produced at all. They take over the entire population very quickly and cause a danger for premature convergence. The spin roulette scheme also has low selection pressure when fitness values are near each other. As a result of these, the average fitness of the new population will usually exceed that of the old, since the population will have been shorn of its weaker members [13, 14]. Avoiding this, tournament selection is used. In tournament selection a number tour of individuals is chosen randomly from the population and compared. The best individual from this group becomes a member of new population [15]. This process is repeated till all members of new population are formed. This selection scheme produces more uniform members among old population while protecting the stronger strings.

3.1.2. Crossover

The reproduction process only discards the poorer strings. Changing of the strings is performed by the crossover operator. It provides exchanging information between randomly selected two parents. A cut-off point is determined and the strings are then cut to produce two head segments and two tail segments. The tail segments are swapped over to generate two new full-length offsprings. After mating process the new strings share some characteristics of each parent. The strings undergo crossover are decided by a crossover probability, which is typically between 0.6 and 1.0.

The traditional method of crossover may produce infeasible offspring if the gene code is decimal numbers. The method can result some code numbers appearing twice in offspring after crossover. In order not to carry same genes for the new chromosome a special process is introduced. First, the right segment of the first parent is kept unchanged to form the right segment of the new chromosome. Second, different codes from the currently formed segment are found and selected from the second parent to form the left segment of the new chromosome. This selection is carried out according to the order in which the codes appear in the second parent. After new chromosome is generated, same process is applied for the second new chromosome by exchanging the roles between parents. This special process maintains the genetic codes from both parents.

The crossover operator discussed above is one point crossover. Since it has restriction when it cannot combine certain features encoded on chromosomes, use of two points crossover operator becomes a solution. After selection of two

parents from the population, two crossover points in a parent are randomly generated based on the string length. The first child is generated by arranging the elements of the selected segment of first parent according to the order in which they appear in the second parent [16]. The remaining elements of first child are held as same as in the first parent. The same procedure is repeated for generation of the second child (see Figure 4).

4 6 1 3 8 10 2 11 12 5 7 9 Parent 1 5 8 7 12 Î 11 3 6 1 t 2 9 10 4 Parent 2 1st Î cut-off point T 2nd cut-off point

Arranging the selected elements 10, 2,11, 12 in the order of Parent 2: 12,

11, 2,10

Arranging the selected elements 3, 6, 1, 2 in the order of Parent 1: 6,1, 3,

2

4 6 1 3 8 12 11 2 10 5 7 9 Child 1

5 8 7 12 11 6 1 3 2 9 10 4 Child 2

Figure 4. Two points crossover with decimal numbers.

3.1.3. Mutation

After crossover operator, a mutation operator is employed to new generation. Mutation operator raises the disparity between the generations and protects the loss of some potential valuable strings. It is implemented by a pairwise swap between the positions of randomly selected two gene codes in the chromosome. Mutation is done with a small mutation probability typically between 0.01 and 0.05.

3.2. Simulated annealing

The concept of SA method was proposed by Metropolis et. al [17] and it was first applied to an optimization problem by Kirkpatrick et. al [18]. SA is based on the analogy between the statistical mechanics and solving of large optimization problems. The behavior of solid-state physics is with many degrees of freedom in thermal equilibrium at a finite temperature. Annealing is a thermodynamic process by which solids are heated and cooled gradually to a crystalline state with minimum energy. Firstly, a solid is heated up to a maximum value at which all molecules of the solid randomly arrange

themselves in the liquid phase. Then, if the liquid is cooled slowly enough, the atoms are able to form themselves in the low energy ground state of a corresponding lattice. This is the optimal solution in a mathematical optimization problem. However, if the liquid is cooled fast, the solid will not attain the ground state, rather it will be a higher energy state that corresponds to a local or suboptimal solution in the mathematical sense.

Simulation of the physical cooling for a fixed value of the temperature T can be done with the Metropolis algorithm [17]. It generates sequences of states of the solid that are equal to the solutions in optimization problem. Given a current state si of the solid with energy Esi: which is the cost of solution, a

small perturbation is applied to the solid and the next state Sj is generated with energy Esj. If the difference in energy, AE = Esj-Esi, is negative, the sate Sj is

accepted as the current state and the process is carried out with the new state. Otherwise, the state is still acceptable providing a probability given by Eq. (2).

prob = e

HBF

where kB is Boltzmann constant which is used to cope with different materials

in thermodynamics process. This allows escaping from local minima by accepting higher energy state. At low temperatures there is still a chance of a system being in a high energy state. The acceptance probability is proportional to the temperature T, which is decreased gradually as the algorithm continues. This is repeated until the system freezes into a steady state. The outline of the algorithm is shown in Figure 5.

1. Select a starting solution si, si e S

2. Set an initial temperature, T > 0 3. Set a temperature reduction factor, X 4. REPEAT

Randomly select Sj e N(si)

AE = E(Sj) - E(si)

IF AE < 0; THEN s. = s_i f

END

ELSE IF exp((-AE/kBT) > Random(0, 1)) THEN s.

END

Set T = X (T) UNTIL Stopping criterion

3.2.1. Configuration

Representing of a possible solution starts with construction of a random solution which forms initial starting configuration [19]. A configuration in the problem is a permutation of the number of assembly points, P, in a body. In the solution space the system should move one state to another. This is achieved by defining a neighborhood. It is a set of feasible solutions that the system can move the other states from a certain state. For a given solution, si,

the neighborhood structure, N(si), is build by pairwise interchange operator. In the method, a solution is obtained by switching the assembly points in positions i and j. The total neighborhood is composed of \N(Si)\=P(P-1)/2 neighbors (see Figure 6). The neighborhood could be changed as the annealing process continues. To access to almost any neighborhood state, starting temperature must be hot enough.

2 5 1 4 3 Sample Initial configuration 2-5 2-1 2-4 2-3 Swaps 5-1 5-4 5-3 1-4 1-3 4-3 Total neighborhood =P(P-1)/2 10 Figure 6. Pairwise exchange operation.

3.2.2. Energy and temperature

The quality of a solution is measured by the energy of certain state. The positions of the assembly points one to another determine the energy of the state in the problem. This is simply the total distance between the assembly points to be visited. It is calculated at every iteration of the algorithm. Once the starting temperature set high enough, the system allows enough iteration at each temperature state. This is for stabilizing the system at relative temperature. The annealing process requires the temperature decrement. It is done using temperature control parameter, which is typically in the range of 0.80-0.99. The algorithm is greatly affected by the initial temperature and temperature reduction rate by accepting or rejecting certain moves.

3.3. Tabu search

TS originally proposed by Glover [20] as an optimization procedure. Being an iterative improvement procedure, TS requires an initial solution and endeavors to determine a better solution in means of escaping local optima based on a short-term memory of recent solutions [21]. The recent solutions are stored in a tabu list, which is a historical record of investigated feasible solutions, and it is utilized to explore new regions that it might not otherwise explore. Thus, tabu restrictions have a role in constraining the search space. Moreover, TS allows backtracking to previous solutions, which may lead to better solutions via a different direction [22]. This is done by using aspiration criterion. Aspiration criterion has a role in guiding the search process. In Figure 7, the outline of the tabu search algorithm is given.

1. Initialisation

Select a starting point st e S

Record the best current solution as sbest = si

Set objective function f(sbest) = f(si)

2. Construct a move to create sj, sj e N(Si), (Si, neighborhood of si)

Evaluate all sj IFf(sj) < Cost_best_move

IF move NOT Tabu OR Aspiration Level > Threshold Cost_best_move = f(sj)

Best_move = sj si = Best_move

END END

Update Tabu_List AND Aspiration Level Update Best

3. IF Stopping criterion = False THEN Go to 2 4. ELSE Return to Best.

Figure 7. Outline of the tabu search algorithm.

For a given combinatorial optimization problem, solution space S constitutes a set of feasible solutions. Each solution si in S hosts a subset of S, Si,

comprising of feasible solutions that may be attained from si in a transition. The feasible solutions reached from si, Si, are considered to neighborhood of si.

TS starts from a randomly generated feasible solution, si, and moves to solution sj, which is the best solution in the neighborhood Si, and then new

iterative procedure takes the best move in the neighborhood as the new current move, even though it might be worse than the current move. A stopping criterion based on a maximum number of moves or a time limit finalizes the procedure with recording the overall minimum solution and the cost.

In addition to the basic version of a TS algorithm that uses short-term memory, there is also a second kind of memory, which is long-term memory that is designed to diversify the search by penalizing the moves, which are most frequently played.

3.3.1. Tabu list

Some moves are not allowed since they are recorded as forbidden in a tabu list. Items placed on such a list remain there for a fixed number of iterations are managed by the size of list. Accepted solutions form the list should be considered as tabu for seven moves, as suggested by Glover [20]. After these successive moves, tabu status a move is removed on the basis of first-in-first-out, and then it becomes again allowable. Building of a tabu list prevents cycling of a move. The length of tabu list affects the final result. Long list may result lower overall minimum cost while it requires more computational time. Again, too short list causes high probability of cycling.

3.3.2. Aspiration

Aspiration is a heuristic rule to release temporarily tabu status of a solution based on cost values. Aspiration helps facilitate flexibility of the procedure while escaping local optima and preventing the search being trapped within a subset of a solution. This rule basically states that at all times during the search, associated with each cost value, c; there is also an aspiration value which is the current lowest cost of all solutions reached from solutions with cost c. In other words, if the cost of a tabu solution is less than the aspiration value related to the cost of the current solution, then the tabu status of the tabu solution is temporarily neglected. Tabu solution is not removed from the tabu list but just ignored and the move is performed [23]. As move is made, the aspiration value is updated providing that the cost of the solution to which the search has advanced is less than the current aspiration value.

In a concise manner, the algorithm requires an initial solution, si, which is a permutation of the number of assembly points to be visited, and the associated cost is c(si). The neighborhoods of each successive solution can be defined as

same as in SA. Among the neighborhood, a solution sng with minimum cost cng

(Cost_best_move) is accepted providing that it is not in the tabu list or its

aspiration value is less than the current solution of sng (Figure 7). It is

should be considered as the standard order of the move. As the algorithm progressed, the current overall minimum solution smin (Best_move) and the

current overall minimum cost cmin (Best) are recorded at each step. After the

specified number of moves the procedure gives the overall minimum solution and its cost value.

4. COMPUTATIONAL RESULTS AND COMPARISON OF

ALGORITHMS

The tests have been made on a P4 2.6 GHz PC and the programming language used has been C. The input data of the coordinates for each assembly point on the base part has been obtained from the computer database. For the GA, the population size was 200 per generation and the number of generations was 100. The probabilities of crossover and mutation were chosen 0.9 and 0.05 respectively and two points crossover operator was employed. The tournament size was set to 5. For the SA, the temperature was reduced at 100 iterations with temperature control value of 0.95. The initial temperature value was chosen 100. For each temperature level, the number of configurations to be produced was decided randomly. In fact, these state repetitions at each temperature can be reduced. This is done for avoiding saturation as it occurs if the cost of function does not change at that temperature. The major drawback of SA is that finding high quality solutions demands large computational effort. The TS algorithm employed short-term memory and aspiration criterion. As suggested, the length of tabu list was chosen seven in the problem and the maximum number of moves was defined as 1500. Implementing of these features meets the aim of the problem, hence other features of the TS algorithm, which are intensification and diversification strategies, have not been implemented.

Each algorithm has been run several times to compute the shortest path of 15 and 30 assembly point problems. A sequencing problem with 15 nodes has roughly 1012 feasible solutions. In other words for a 30 nodes problem makes

30! ^ 1032 possible combinations. These imply that because of the search space

is quite big it would take very long time to find optimal solution with classical methods. On the other hand, the algorithms used have converged to the best objective function values quickly shown in Figure 8 and Figure 9. It can be seen that GA, TS and SA can converge to same optimal solution. Although GA is the most effective algorithm regarding having the quickest converge time. For the case of 15 assembly points, the best objective function values have been obtained at the iteration of 9 for GA, iteration 17 for TS and the iteration 19 for SA algorithms (Figure 8). For the case of 30 assembly points, the best function convergences for GA, TS and SA were at the iterations of 15, 81 and 77 respectively (Figure 9). The algorithms were stopped when there was no

changes at the function values after a number of iterations. The GA provided to be faster than SA and TS algorithms. The rest algorithms were converged almost at the same iterations. It can be seen that when the size of the problem increases the time required for the convergence also increases.

I t e r a t i o n s

Figure 8. Best function convergences for 15 assembly points.

Iterations

5. CONCLUSIONS

In order to practice automation of assembly sequences as a part of CIM process, a prototype of gantry-type robot was designed and produced. A self-reasoning system about robot movements was integrated with the prototype. The system acquires the input data, coordinates of the parts to be assembled, from computer database. Then, three metaheuristic techniques, GA, SA and TS, engage to find the shortest path of assembly sequence. After comparison of the results the system takes the optimum assembly configuration to pass automatically to the control unit for carrying out assembly tasks.

The problem of optimum sequencing in an assembly automation process is outlined. The problem states that a number of parts should be assembled such a sequence that a shortest assembly route can be achieved. The parts to be assembled are assumed having same assembly characteristics so that the robot gripper can handle those parts at same condition. Assembly sequence is randomly generated and each part is assembled at once.

The proposed system obtains quick solutions depending on the number of parts assembled in assembly process. The algorithms are shown to be effective regarding the results obtained. The automation of assembly sequences that are stressed as the aim of study has been realized despite the performance differences among the algorithms. In order to increase precision of SA and TS algorithms, some further work may be needed. However, the model developed can be employed as a self-reasoning system for facilitating automation of assembly in a CIM environment.

References

[1] Baartman, J. P.; Automation of Assembly Automations on Parts, PhD Thesis, Delft University of Technology, The Netherlands (1995).

[2] Asfahl, C. R.; Robots and Manufacturing Automation, 2nd ed., John-Willey & Sons Inc., New York (1992).

[3] Bhaskare, A. and Newton, R.; Assembly sequence design methods, Proceedings of the International Conference on Assembly, Adelaide (1993), 109-115. [4] Boneschanscher, N.; Plan Generation for Flexible Assembly Cells, PhD Thesis,

Delft University of Technology, The Netherlands (1993).

[5] Hormann, A. and Rembold, U.; Development of an advanced robot for autonomous assembly, Proceedings of IEEE Int. Conf. on Robotics and Automation, Sacramento, California (1991) 2452-7.

[6] Rohatynski, R., Kielec, R., Sasiadek, M.; A computer assisted assembly sequence planning in mechanical design, AI-METH Artificial Intelligence Methods. Gliwice, Poland (2003), 112-3.

[7] Rehg, J. A. and Kraebber, H. W.; Computer Integrated Manufacturing, 2nd ed., Prentice-Hall, New Jersey (2001).

[8] Gologlu, C.; Intelligent short path planning for assembly sequence automation in CIM, Proceedings of IMS'2004 4th International Symposium on Intelligent Manufacturing Systems. Sakarya, Turkey (2004), 290-7.

[9] Costamagna, E., Fanni, A., Giacinto, G.; A tabu search algorithm for the optimization of telecommunication networks, European Journal of Operational Research, 106 (1998), 357-372.

[10] Mergenthaler, W., Stadler, W., Wilbertz, H., Zimmer, N.; Optimizing automotive manufacturing sequences using simulated annealing and genetic algorithms, Control Engineering Practice, 3(4) (1995), 569-573.

[11] Pham, D. T., Karaboga, D.; Intelligent Optimisation Techniques, Springer-Verlag, London (2000).

[12] Holland, J. H.; Adaptation in Natural and Artificial Systems, University of Michigan Press, Michigan (1975).

[13] Eiben, A. E., Smith, J. E.; Introduction to Evolutionary Computing, Springer, The Netherlands (2003).

[14] Cartwright, H. M.; The genetic algorithm in science, Pesticide Science, 45 (1995), 171-178.

[15] Reeves, C. R.; Modern heuristic techniques, In: Rayward-Smith, V. J., Osman, I. H., Reeves, C. R., Smith, G. D., editors, Modern heuristic search methods, John-Willey & Sons, England (1996), 1-25.

[16] Bhaskara, Reddy S. V., Shunmugam, M. S., Narendran, T. T.; Operation sequencing in CAPP using genetic algorithms, International Journal of Production Search, 37(5) (1999), 1063-1074.

[17] Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A., Teller, E.; Equations of state calculations by fast computing machine, Journal of Chemical Physics, 21 (1953), 1087-1091.

[18] Kirkpatrick, S., Gelatt, Jr. C. D., Vecchi, M. P.; Optimization by simulated annealing, Science, 220 (4598) (1983), 671-680.

[19] Khan, Z., Prasad, B., Singh, T.; Machining condition optimization by genetic algorithms and simulated annealing, Computers and Operational Research, 24 (7) (1997), 647-657.

[20] Glover, F.; Heuristics for integer programming using surrogate constraints, Decision Sciences, 8 (1977) 156-166.

[21] Thesen, A.; Design and evaluation of tabu search algorithms for multiprocessor scheduling, Journal of Heuristics, 4 (1998) 141-160.

[22] Winker, P., Gilli, M.; Applications of optimization heuristics to estimation and modelling problems, Computational Statistics and Data Analysis, 47 (2) (2004), 211-223.

[23] Bland, J. A., Dawson, G. P.; Tabu search and design optimization, Computer-Aided Design, 23 (3) (1991), 195-201.