This is an HTML rendering of a working paper draft that led to a publication. The publication should always be cited in preference to this draft using the following reference:

Citation(s): 50 (Google Scholar), 33 (selected).

This document is also available in PDF format.

The document's metadata is available in BibTeX format.

Find the publication on Google Scholar

This material is presented to ensure timely dissemination of scholarly and technical work. Copyright and all rights therein are retained by authors or by other copyright holders. All persons copying this information are expected to adhere to the terms and constraints invoked by each author's copyright. In most cases, these works may not be reposted without the explicit permission of the copyright holder.

Diomidis Spinellis Publications

Stochastic Algorithms for Buffer Allocation in Reliable Production Lines

Diomidis D. Spinellis Chrissoleon T. Papadopoulos
Department of Information and Communication Systems Department of Business Administration
GR-832 00 Karlovassi GR-821 00 Chios
University of the Aegean University of the Aegean
Greece Greece


The allocation of buffers between workstations is a major optimization problem faced by manufacturing systems designers. It entails the determination of optimal buffer allocation plans in production lines with the objective of maximizing their throughput. We present and compare two stochastic approaches for solving the buffer allocation problem in large reliable production lines. The allocation plan is calculated subject to a given amount of total buffer slots using simulated annealing and genetic algorithms. The throughput is calculated utilizing a decomposition method.


Simulated annealing, genetic algorithms, buffer allocation

AMS Classification:

90C15 Stochastic programming

1  Introduction

The allocation of buffers between workstations is a major optimization problem faced by manufacturing systems designers. It has to do with devising an allocation plan for distributing a certain amount of buffer space among the intermediate buffers of a production line. It is a very complex task that must account for the random fluctuations in mean production rates of the individual workstations of the lines. To solve this problem there is a need for two different tools. The first is a tool that calculates the performance measure of the line which has to be optimized (e.g., the throughput or the mean work-in-process). This may be an evaluative method such as simulation, a decomposition method [1,2], or a traditional Markovian state model in conjunction with an exact numerical algorithm [3]. The second tool is a search (generative) method that tries to determine an optimal or near optimal value for the decision variables, which in our case are the buffer capacities of the intermediate buffer locations in the line. Examples of such methods are the classical search methods such as the well-known Hooke-Jeeves method, various heuristic methods, knowledge based methods, genetic algorithms, and simulated annealing.

Evaluative and generative (optimization) models can be combined in a `closed loop' configuration by using feedback from an evaluative model to modify the decision taken by the generative model. In such a configuration the evaluative model is used to obtain the value of the objective function for a set of inputs. The value of the objective function is then communicated to the generative model which uses it as an objective criterion in its search for an optimal solution. In the rest of this paper we will use the formalism S(G, E) to describe a closed loop system using the generative method G and the evaluative method E. The generative models that will be used in this paper are:

complete enumeration,
reduced enumeration,
genetic algorithms, and
simulated annealing.

Furthermore, two evaluative models will be used:

the exact numerical algorithm described in the work by Heavey at al. [3], and
the decomposition algorithm numbered as A3 in the work by Dallery and Frein [2].

An overview of the existing research in the area of evaluative and generative models of manufacturing systems can be found in two review papers [4,5] and a number of books [6,7,8,9,10,11].

Several researchers have studied the problem of optimizing buffer allocation to maximize the efficiency of a reliable production line [12,13,14]. Their results are based on comprehensive studies to characterize the optimal buffer allocation pattern. Authors have provided extensive numerical results for balanced lines with up to 6 stations and limited results for lines with up to 9 stations. However, few methods can handle this problem for large production lines, in a computationally efficient way. In this paper we compare two stochastic approaches suitable for large production lines, one based on genetic algorithms, and one based on simulated annealing. Details on how these methods can be applied to the problem can be found in the work by Bulgak et al. [15] which describes the application of genetic algorithms for the buffer allocation in asynchronous assembly systems and in our work [16] which describes a corresponding approach using simulated annealing. The implementation of both approaches in this paper works in close cooperation with a decomposition method [2].

Simulated annealing is an adaptation of the simulation of physical thermodynamic annealing principles [17] to the combinatorial optimization problems [18,19]. Similar to genetic algorithms and tabu search techniques [20] it follows the ``local improvement'' paradigm for harnessing the exponential complexity of the solution space. The algorithm is based on randomization techniques. An overview of algorithms based on such techniques can be found in the survey by Gupta et al. [21]. A complete presentation of the method and its applications is described by Van Laarhoven and Aarts [22] while a number of works present accessible algorithms for its implementation [23,24]. As a tool for operational research simulated annealing is presented by Eglese [25], while Koulamas et al. [26] provide a complete survey of simulated annealing applications to operations research problems.

Genetic algorithms [27,28,29] are global optimization techniques that avoid many of the shortcomings exhibited by local search techniques on difficult search spaces, such as the buffer allocation problem. Goldberg [30] describes a number of diverse genetic algorithm applications, while Karr [31] presents their use for modelling, design, and process control. Finally, Tompkins and Azadivar [32] use genetic algorithms for optimizing simulated systems.

This paper is organized as follows: Section 2 states the problem and the assumptions of the model and Section 3 describes the evaluation methodology and associated implementation decisions. In Section 4, we compare the numerical results obtained from the algorithms. Finally, Section 5 concludes the paper and suggests some future research directions.

2  The Buffer Allocation Problem

In asynchronous production lines, each part enters the system from the first station, passes in order from all stations and the intermediate buffer locations, and exits the line from the last station. The flow of the parts works as follows: in case a station has completed its processing and the next buffer has space available, the processed part is passed on. Then, the station starts processing a new part that is taken from its input buffer. In case the buffer has no parts, the station remains empty until a new part is placed in the buffer. This type of line is subject to manufacturing blocking (or blocking after service) and starving.

2.1  Assumptions of the Model

The model operates under the assumption that the first station is never starved and the last station is never blocked. The processing (service) times at each station are assumed to be independent random variables following the exponential distribution, with mean service rates, mi, i = 1,2,, K. In our model, the stations of the line are assumed to be perfectly reliable, that is, breakdowns are not allowed.

The exponentiality of the processing times as well as the absolute reliability of the line's workstations are rather unrealistic assumptions. However, the service completion times can be exponential or can be approximated by an exponential distribution. The variability in completion times may be attributed to failures and repairs which implicitly exist in the problem at hand. Following this view, the proposed model may be applied to any unreliable production line under the exponentiality assumptions for the service completion times.

Figure 1 depicts a K-station line that has K-1 intermediate locations for buffers, labeled B2, B3, , BK.

A K-station production line with K-1 intermediate buffers
Figure 1: A K-station production line with K-1 intermediate buffers

The basic performance measures in the analysis of production lines are the throughput (or mean production rate) and the average work-in-process (WIP) or equivalently the average production (sojourn) time.

The object of the present work is the buffering of asynchronous, reliable production lines with the assumptions given above. The objective is the maximization of the line's throughput, subject to a given total buffer space.

2.2  The Buffer Allocation Problem

In mathematical terms, the buffer allocation problem can be stated as follows:

Find B = (B2,B3,,BK) so as to



subject to:


i = 2 
Bi = N
Bi 0
Bi   integer  (i = 2,3,,K)


The number of feasible allocations of N buffer slots among the K-1 intermediate buffer locations increases dramatically with N and K and is given by the formula:

N + K - 2
K - 2

= (N+1)(N+2)(N+K-2)

For this reason exhaustive search techniques are not practical for determining optimal configurations of production lines with a large number of stations or buffers.

3  Evaluation Methodology

We have evaluated different approaches for solving the optimal buffer allocation problem for large production lines, by performing the following steps:

We utilized the decomposition method [2] as an evaluative tool to determine the throughput of the lines. The algorithm computes approximately the throughput for any K-station line with finite intermediate buffers and exponentially distributed processing times.

To find the buffer allocation that maximizes the throughput of the line, we utilized two stochastic methods, simulated annealing and genetic algorithms, specifically adapted for solving this problem.

In order to evaluate the applicability of the stochastic methods to the buffer allocation problem using comparable architectures we designed and implemented a system to calculate the optimum buffer configuration for a given reliable production line using a variety of algorithms [33]. The system takes as input:

Based on the above input, the system calculates the buffer allocations B = (B2,B3,,BK) for the maximal line throughput. Furthermore, the system is instrumented to provide as part of the solution the throughput of the suggested configuration, as well as the number of different configurations that were tried. The line throughput is used to evaluate the quality of the suggested configuration when compared with the throughput calculated by other methods. The number of different configurations tried, is used as an objective performance criterion, because the configuration evaluation step is the dominant execution time factor and the basic building block of all optimization methods. In addition, a special system configuration allows the creation of a file containing step-by-step snapshots of the algorithm progress. After obtaining the test results we wrote a number of scripts in the Perl programming language [34] that utilized the snapshot file to visualize and animate the dynamic behavior of the algorithms.

We ran a number of tests on both balanced and unbalanced lines and compared the stochastic method results against each other and against the results obtained by other methods. For short lines and limited buffer space a complete enumeration of all configurations provided an accurate measurement base to verify the stochastic algorithm results. For larger configurations we used a reduced enumeration in order to provide the comparative measure.

3.1  The Reduced Enumeration Method

Reduced enumeration is based on the experimental observation that the absolute difference of the respective elements of the optimal buffer allocation (OBA) vectors with N and N+1 buffer slots is less than or equal to 1:

|BiN+1 - BiN| 1,   "i: 2 i K.

In this way, we have been able to derive the OBA by induction for any number N of buffer slots that are to be allocated among the K-1 buffer locations of the line. The reduction works as follows: when N* and K are given one needs to determine all the OBA vectors for N = 1,2,,N* and then for N = N*+1 by searching only the values of BiN-1, BiN and BiN+1. Furthermore, this reduction starts after a number of total buffer slots N. To quantify the reduction, by applying the improved enumeration it has been experimentally observed that the number of iterations were reduced by at least 60% for short lines. This reduction accounts for well over 90% for large production lines (with more than 12 stations). Recall that the number of feasible allocations of N buffer slots among the K-1 intermediate buffer locations increases dramatically with N and K and is given by Formula (3).

3.2  Simulated Annealing

Simulated annealing is an optimization method suitable for combinatorial minimization problems. Such problems exhibit a discrete, factorially large, configuration space. In common with all paradigms based on ``local improvements'' the simulated annealing method starts with a non-optimal initial configuration (which may be chosen at random) and works on improving it by selecting a new configuration using a suitable mechanism (at random in the simulated annealing case) and calculating the corresponding cost differential (D[`O]K). If the cost is reduced, then the new configuration is accepted and the process repeats until a termination criterion is satisfied. Unfortunately, such methods can become ``trapped'' in a local optimum that is far from the global optimum. Simulated annealing avoids this problem by allowing ``uphill'' moves based on a model of the annealing process in the physical world.

Our implementation of the simulated annealing algorithm for distributing N buffer space in a K-station line [16] follows the following steps:

  1. [Set initial line configuration.] Set Bi N/K , set BK/2 BK/2 + N - i = 2K N/K .
  2. [Set initial temperature T0.] Set T 0.5.
  3. [Initialize step and success count.] Set S 0, set U 0.
  4. [Create new line with a random redistribution of buffer space.] Move Rn space from a source buffer BRs to a destination buffer BRd: set B B, set Rs rand [2 K+1), set Rd rand [2 K+1), set Rn rand [0 BRs+1), set BRs BRs - Rn, set BRd BRd + Rn.
  5. [Calculate energy differential.] Set DE [`O]K(B) - [`O]K(B).
  6. [Decide acceptance of new configuration.] Accept all new configurations that are more efficient and, following the Boltzmann probability distribution, some that are less efficient: if DE < 0 or exp([(-DE)/ T]) > rand (0 1), set B B, set U U + 1.
  7. [Repeat for current temperature]. Set S S + 1. If S < maximum number of steps, go to step 4.
  8. [Lower the annealing temperature.] Set T c T (0 < c < 1).
  9. [Check if progress has been made.] If U > 0, go to step 3; otherwise the algorithm terminates.

3.3  Genetic Algorithms

Genetic algorithms are also global optimization techniques that avoid many of the shortcomings exhibited by local search techniques on difficult search spaces. They rely on modeling the problem as a population of organisms. Every organism represents a possible valid solution to the problem. Organisms are composed of alleles representing parts of a given solution. Standard genetic recombination operators are used to create new organisms out of existing ones by combining alleles of the existing organisms. In addition, mutations can randomly change the composition of existing organisms. Typically, the algorithm evaluates all the organisms of the population and creates new organisms by combining existing ones based on their fitness. This procedure is repeated until the variance of the population reaches a predefined minimum value.

An important characteristic of our implementation of the genetic algorithm concerns the representation of the solution. A good representation should ensure that the application of standard crossover recombination operators (where a new organism is composed from parts of two existing ones) would result in a valid new representation. Representing the line configuration as a vector B of buffers allocated across the line is not such a representation since given two buffer configurations (B1, B2) and recombining them as a new buffer B at point c so that B0 c B1, 0 c and Bc + 1 K B2, c + 1 K will not guarantee that i = 1KBi = N i.e. that the resulting line configuration will be composed of N buffers. For this reason we devised an alternative, position-based, representation using a vector P of length equal to the number of buffers N. Every element of P can take values 0 K representing the position of the given buffer slot within the production line. The two representations are equivalent; the vector P can be mapped to B as follows:

Bi = N

j = 0 

if i = j

The position-based representation will generate valid buffer configurations using standard genetic crossover and mutation operators. Using this representation, the genetic algorithm we implemented for distributing N buffer space in a K-station line is described in the following steps:

  1. [Initialize a population of size S.] Set P0 S,0 N rand[0 K-1).
  2. [Evaluate population members creating throughput vector T.] For i 0 S: set Ti [`O]K(Pi).
  3. [Create roulette selection probability vector R . ] Set Ri j = 0i (Tj / k = 0STk).
  4. [Create new population using crossovers from the previous population.] For i 0 S: if rand [0 1) < crossover rate, set c rand [0 S), set Pi, 0 c PRr, 0 c, set Pi, c + 1 N PRr, c + 1 N; otherwise set Pi PRr by selecting each r using the roulette selection probability vector so that Rr rand [0 1) < Rr+1.
  5. [Introduce mutations.] For i 0 S: for j 0 N: if rand [0 1) < mutation rate, set Pi,j rand[0 K-1).

  6. [Keep fittest organism for elitist selection strategy.] Select f so that Tf T0 S, set Prand [0 S) Pf.
  7. [Make new population the current population.] Set P P.
  8. [Loop based on the population's variance.] If i = 0P | Tf - Ti | > maximum variance go to step 2; otherwise the algorithm terminates with the optimal line setup in Pf.

The implementation of genetic algorithms can be tuned using a number of different parameters. In our implementation we used the parameters that Grefenstette [35] derived using meta-search techniques namely:

The random floating point numbers 0 < R < 1 used for selecting energy differentials based on the annealing temperature R < exp([(-DE)/ T]), the crossover points, the mutation rates, and the selection of organisms are produced using the subtractive method algorithm [36]. Finally, the evaluative function that we used for calculating DE is based on the decomposition method [2].

4  Method comparison

Figure 2: Stochastic method operation comparison

Before detailing the comparative results of our examination, it is interesting to visualize the operation of the two stochastic methods. Figure 2 depicts the runtime behavior of the two methods. Each point on the two scatter charts represents a given production line throughput value at a specific step of the algorithm. Both charts depict the calculation of the placement of 30 buffers in a balanced line of 10 stations. The simulated annealing algorithm optimizes a single solution in the specific example in 45.000 iterations. The solution's throughput value oscillates as both better and worse solutions are randomly selected at each iteration step. As can be seen on the chart, the oscillation width decreases following the algorithm's exponential cooling schedule and converges towards the optimal value.

The genetic algorithm is based on the implicit parallelism of the solutions represented by the initial population depicted on the chart's z-axis. Thus, in the specific example, it terminates with an optimal configuration after 250 generations. As the chart demonstrates the search starts with a wide spectrum of different solutions which are evaluated and evolve in parallel with non-optimal solutions gradually becoming extinct. Mutations and recombinations regenerate suboptimal solutions, but, due to the ``survival of the fittest'' organism selection strategy, their survival does not last for long.

Figure 3: Computed throughput of simulated annealing S(SA, Deco) and genetic algorithms S(GA, Deco) compared with complete enumerations using the exact S(CE, Exact) and the decomposition evaluative methods S(CE, Deco) for 9 stations (left); compared with reduced enumeration S(RE, Deco) for 15 stations (right).

Our first comparison experiment concerned the algorithm operation on balanced lines for cases where exact solutions were known. In Figure 3 we present the optimum throughput configurations for balanced lines found using the stochastic methods against the throughput found using complete (for 9 stations) and reduced enumeration techniques. It is apparent that the stochastic algorithm results are almost identical and follow closely the results obtained by the other methods. Both methods are subject to the reduced evaluative accuracy of the decomposition method compared to the Markovian model.

Figure 4: Simulated annealing S(SA, Deco) and genetic algorithms S(GA, Deco) with decomposition evaluation versus complete enumerated Markovian S(CE, Exact) throughputs for unbalanced lines with 4-6 stations.

In addition to the balanced line evaluation, we compared the stochastic methods against unbalanced line enumeration using the Markovian evaluative procedure for a variety of line sizes, service time configurations, and available buffer space. The results are summarized in Figure 4. It is apparent, that the stochastic method configurations - although identical to each other - are not always optimal for limited available buffer space; however, they quickly converge towards the optimal configurations as buffer space increases. This difference can be accounted by the use of the fast decomposition evaluative procedure used in the stochastic algorithm implementation yielding approximate results against the use of the Markovian evaluative procedure for the enumeration method yielding exact results.

Figure 5: Performance of simulated annealing S(SA, Deco) and genetic algorithms S(GA, Deco) compared with complete S(CE, Deco) and reduced S(RE, Deco) enumerations for 9 stations and 15 stations. Note the log10 scale on the ordinate axis.

Our goal for using stochastic methods is to optimize large production line problems where the cost of other methods is prohibitively expensive. As an example the reduced enumeration method when run on a 15 station line with a buffer capacity of 30 units took more than 10 hours to complete on a 100MHz Pentium processor. As shown in Figure 5 the cost of the stochastic methods is higher than the cost of the full and reduced enumeration methods for small lines and buffer allocations. However, it quickly becomes competitive as the number of stations and the available buffer size increase. In addition, the performance of the genetic algorithm implementation is approximately an order of magnitude better than the simulated annealing implementation. Notice that - in contrast to the deterministic methods - the stochastic method cost does not increase together with the available buffer space and that it increases only linearly with the number of stations.

Figure 6: Performance and accuracy of simulated annealing S(SA, Deco) compared with genetic algorithms S(GA, Deco) for large production lines. Note the log10 scale on the ordinate axis.

Finally, Figure 6 depicts the comparative performance and calculated throughput for the two stochastic methods when optimizing lines of up to 400 stations and 1200 buffers. The genetic algorithm implementation producing solutions with only 2.000.000 evaluations even for 400 station lines is clearly the performance winner. However, as depicted on the right hand chart, the throughput of the line configuration found by the genetic algorithm is consistently lower than the throughput of the line found by the simulated annealing method. The results we obtained could not be independently verified, because no other numerical results for the buffer allocation problem in large production lines can be found in the open literature.

5  Conclusions and Future Directions

The results obtained applying stochastic methods to the reliable line near-optimal buffer allocation problem are interesting. The performance and the accuracy of the methods, although inferior for optimizing small lines with limited buffer space, indicate clearly that they become the methods of choice as the problem size increases. Both methods can be used for optimizing large line configurations with simulated annealing producing more optimal configurations and the genetic algorithm approach leading in performance. This indicates that the two methods can be used in complimentary fashion. Real-time applications can utilize genetic algorithms for the swift recalculation of optimal configurations, while batch-oriented calculations can utilize simulated annealing for obtaining an optimal configuration.

Further investigation is needed in order to fully evaluate the potential of the two methods. The failure, in large production lines, of the genetic algorithm method to locate the optimal configuration found by the simulated annealing method is intriguing. It would be interesting to carefully examine the ``endgames'' of the two methods and find if and how the genetic algorithm implementation can be tweaked to evolve towards more optimal configurations. A dynamic re-adjustment of the algorithm's parameters (population size, crossover rate, mutation rate, etc.) forms one such possibility.

The annealing schedule and the genetic algorithm parameters that we used can clearly be optimized potentially increasing both methods' accuracy and performance. The use of heuristics in setting up the initial buffer configuration can decrease the number of steps needed for reaching the optimal. The differing relative strengths of the two stochastic approaches could also be combined in the form of a hybrid algorithm. Such an algorithm could capitalize on the rapid convergence exhibited by the genetic algorithms to quickly arrive at an acceptably efficient solution pruning away dead-ends. It could then pass the quickly derived buffer configuration to a simulated annealing algorithm which would use it as a starting point for obtaining an optimal solution.


S. B. Gershwin. An Efficient Decomposition Method for the Approximate Evaluation of Tandem Queues with Finite Storage Space and Blocking. Oper. Res., 35(2):291-305, 1987.

Y. Dallery and Y. Frein. On Decomposition Methods for Tandem Queueing Networks with Blocking. Oper. Res., 41(2):386-399, 1993.

C. Heavey, H. T. Papadopoulos, and J. Browne. The Throughput Rate of Multistation Unreliable Production Lines. European J. Oper. Res., 68:69-89, 1993.

Y. Dallery and S. B. Gershwin. Manufacturing Flow Line Systems: A Review of Models and Analytical Results. Queueing Systems: Theory and Applications, 12:3-94, 1992.

H. T. Papadopoulos and C. Heavey. Queueing Theory in Manufacturing Systems Analysis and Design: A Classification of Models for Production and Transfer Lines. European J. Oper. Res., 92:1-27, 1996.

H. T. Papadopoulos, C. Heavey, and J. Browne. Queueing Theory in Manufacturing Systems Analysis and Design. Chapman and Hall, London, 1993.

R. G. Askin and C. R. Standridge. Modeling and Analysis of Manufacturing Systems. John Wiley, New York, 1993.

J. A. Buzacott and J. G. Shanthikumar. Stochastic Models of Manufacturing Systems. Prentice Hall, New Jersey, 1993.

S. B. Gershwin. Manufacturing Systems Engineering. Prentice Hall, New Jersey, 1994.

H. Perros. Queueing Networks with Blocking. Oxford University Press, 1994.

T. Altiok. Performance Analysis of Manufacturing Systems. Springer-Verlag, New York, 1997.

F. S. Hillier and K. C. So. The Effect of the Coefficient of Variation of Operation Times on the Allocation of Storage Space in Production Line Systems. IIE Trans., 23(2):198-206, 1991.

F. S. Hillier, K. C. So, and R. W. Boling. Notes: Toward Characterizing the Optimal Allocation of Storage Space in Production Line Systems with Variable Processing Times. Management Science, 39(1):126-133, 1993.

S. B. Gershwin and J. E. Schor. Efficient Algorithms for Buffer Space Allocation. In International Workshop on Performance Evaluation and Optimization of Production Lines, pages 217-228, Samos, Greece, May 1997. University of the Aegean, Department of Mathematics.

A. A. Bulgak, P. D. Diwan, and B. Inozu. Buffer Size Optimization in Asynchronous Assembly Systems Using Genetic Algorithms. Computers ind. Engng., 28(2):309-322, 1995.

D. Spinellis and C. T. Papadopoulos. A Simulated Annealing Approach for Buffer Allocation in Reliable Production Lines. Annals of Operations Research, 1999. To appear.

N. Metropolis, A. N. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. Equation of State Calculation by Fast Computing Machines. J. Chem. Phys., 21(6):1087-1092, 1953.

S. Kirkpatrick, C. D. G. Jr., and M. P. Vecchi. Optimization by Simulated Annealing. Science, 220:671-679, 1983.

V. Cerny. Thermodynamical Approach to the Traveling Salesman Problem: an Efficient Simulation Algorithm. J. Optim. Theory Appl., 45:41-51, 1985.

F. Glover. Tabu Search - Part I. ORSA Journal on Computing, 1:190-206, 1990.

R. Gupta, S. A. Smolka, and S. Bhaskar. On Randomization in Sequential and Distributed Algorithms. ACM Comput. Surv., 26(1):7-86, 1994.

P. J. M. Van Laarhoven and E. H. L. Aarts. Simulated Annealing: Theory and Applications. D. Reidel, Dordrecht, The Nethelands, 1987.

A. Corana, M. Marchesi, C. Martini, and S. Ridella. Minimizing Multimodal Functions of Continuous Variables with the ``Simulated Annealing'' Algorithm. ACM Trans. Math. Software, 13(3):262-280, Sept. 1987.

W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipes in C, pages 343-352. Cambridge University Press, 1988.

R. W. Eglese. Simulated Annealing: A Tool for Operational Research. European J. Oper. Res., 46:271-281, 1990.

C. Koulamas, S. R. Antony, and R. Jaen. A Survey of Simulated Annealing Applications to Operations Research Problems. Omega International Journal of Management Science, 22(1):41-56, 1994.

J. H. Holland. Adaptation in Natural and Artificial Systems. University of Michigan Press, Ann Arbor, Michigan, 1975.

D. E. Goldberg. Genetic Algorithms: In Search of Optimization & Machine Learning. Addison-Wesley, 1989.

S. Forrest. Genetic Algorithms. ACM Comput. Surv., 28(1):77-83, Mar. 1996.

D. E. Goldberg. Genetic and Evolutionary Algorithms Come of Age. Commun. ACM, 37(3):113-119, Mar. 1994.

C. L. Karr. Genetic Algorithms for Modelling, Design, and Process Control. In CIKM '93. Proceedings of the Second International Conference on Information and Knowledge Management, pages 233-238. ACM, 1993.

G. Tompkins and F. Azadivar. Genetic Algorithms in Optimizing Simulated Systems. In WSC '95. Proceedings of the 1995 Conference on Winter Simulation, pages 757-762. ACM, 1995.

D. Spinellis and C. T. Papadopoulos. ExPLOre: A Modular Architecture for Production Line Optimisation. In D. K. Despotis and C. Zopounidis, editors, Proceedings of the 5th International Conference of the Decision Sciences Institute, DSI '99, pages 1446-1449, Athens, Greece, July 1999. Decision Sciences Institute.

L. Wall and R. L. Schwartz. Programming Perl. O'Reilly and Associates, Sebastopol, CA, USA, 1990.

J. J. Grefenstette. Optimization of Control Parameters for Genetic Algorithms. IEEE Transactions on Systems, Man, and Cybernetics, 16(1):122-128, 1986.

D. E. Knuth. The Art of Computer Programming, volume 2 / Seminumerical Algorithms, pages 171-173. Addison-Wesley, second edition, 1981.