Keywords: Buffer Allocation, Non-linear, Stochastic, Integer, Network Design, Simulated Annealing

Introduction and literature review

A large amount of research has been devoted to the analysis and design of production lines. A lot of this research concerns the design of manufacturing systems with considerable inherent variability in the processing times at the various stations, a common situation with human operators/assemblers. The literature on the modelling and optimisation of production lines is vast, allowing us to review only the most directly relevant studies. For a systematic classification of the relevant works on the stochastic modelling of these and other types of manufacturing systems (e.g., transfer lines, flexible manufacturing systems (FMS), and flexible assembly systems (FAS)), the interested reader is referred to a review paper by Papadopoulos and Heavey (1996) and some recently published books, such as those by Askin and Standridge (1993), Buzacott and Shanthikumar (1993), Gershwin (1994), Papadopoulos et al. (1993), Viswanadham and Narahari (1992), and Altiok (1997). There are two basic problem classes:

1. the evaluation of production line performance measures such as throughput, mean flow time, workstation mean queue length, and system utilisation, and

2. the optimisation of the decision variables of these lines.

Examples of decision variables that have been considered are:

1. the sizes of the buffers placed between successive workstations of the lines,

2. the number of servers allocated to each workstation, and

3. the amount of workload allocated to each workstation.

The corresponding optimisation problems are named, respectively, (1) the buffer allocation problem, (2) the server allocation problem, and (3) the workload allocation problem, in a production line.

In Papadopoulos et al. (1993) both evaluative and generative (optimisation) models are given for modelling the various types of manufacturing systems. This work falls into the second category. Evaluative and optimisation models can be combined by closing the loop between them; that is, one can use feedback from an evaluative model to modify the decisions taken by the optimisation model.

One of the key questions that the designers face in a serial production line is the buffer allocation problem (BAP), i.e., how much buffer storage to allow and where to place it within the line. This is an important question because buffers can have a great impact on the efficiency of the production line. They compensate for the blocking and the starving of the line’s stations. For this reason, the buffer allocation problem has received a lot more attention than the other two design problems. Buffer storage is expensive due both to its direct cost, and to the increase of the work-in-process (WIP) inventories. In addition, the requirement to limit the buffer storage can also be a result of space limitations in the shop floor. The literature on the BAP is extensive. A systematic classification of the research work in this area is given in Singh and MacGregor Smith (1997) and Papadopoulos et al. (1993). The works are split according to:

The method used to solve the buffer allocation problem

search methods, Altiok and Stidham (1983), Smith and Daskalaki (1988), Seong, Chang, and Hong (1995), F. S. Hillier and So (1991a), F. S. Hillier and So (1991b), dynamic programming methods, Kubat and Sumita (1985), Jafari and Shanthikumar (1989), Yamashita and Altiok (1997), among others, and simulation methods, Conway et al. (1988), Ho, Eyler, and Chien (1979).

The type of production line

balanced/unbalanced, (Powell (1994) presents a literature review according to this scheme), or reliable/unreliable; the majority of papers deal with reliable lines. Only a few algorithms have been developed to calculate the performance measures of unreliable production lines Glassey and Hong (1983), Choong and Gershwin (1987), Gershwin (1989), Heavey, Papadopoulos, and Browne (1993). Hillier and So (1991a) and Seong, Chang, and Hong (1995) have dealt with the buffer allocation problem in unreliable production lines.

The number of machines available at each workstation

Similarly, few researchers have dealt with the buffer allocation problem in production lines with multiple (parallel) machines at each workstation Hillier and So (1995), Magazine and Stecke (1996), Singh and Smith (1997).

Apart from the buffer allocation problem, the other two interesting design problems have also been considered by some researchers, e.g., the work allocation problem Hillier and Boling (1979), Hillier and Boling (1966), Hillier and Boling (1977), Ding and Greenberg (1991), Huang and Weiss (1990), Shanthikumar, Yamazaki, and Sakasegawa (1991), Wan and Wolff (1993), Yamazaki, Sakasegawa, and Shanthikumar (1992), among others, and the server allocation problem Magazine and Stecke (1996), Hillier and So (1989). Hillier and So (1995) studied various combinations of these three design problems. Other references may be found therein Buzacott and Shanthikumar (1992), Buzacott and Shanthikumar (1993), among others.

The present work deals with the same design problems (buffer allocation, server allocation, and workload allocation) but for long production lines with multi-machine stations.

As the problem being investigated is combinatorial in nature, traditional Operations Research techniques are not as practical for obtaining optimal solutions for long production lines. We propose a simulated annealing (SA) approach as the search method in conjunction with the expansion method developed by Kerbache and MacGregor Smith (1987) as the evaluative tool. Simulated annealing is an adaptation of the simulation of physical thermodynamic annealing principles described by Metropolis et al. (1953) to the combinatorial optimisation problems Kirkpatrick, Gelatt Jr., and Vecchi (1983), Cerny (1985). Similar to genetic algorithms Holland (1975), Goldberg (1989) and tabu search techniques Glover (1990) it follows the ‘local improvement’ paradigm for harnessing the exponential complexity of the solution space.

The algorithm is based on randomisation techniques. An overview of algorithms based on such techniques can be found in the survey by Gupta et al. (1994). A complete presentation of the method and its applications is described by Van Laarhoven and Aarts (1987) and accessible algorithms for its implementation are presented by Corana et al. (1987) and Press et al. (1988). A critical evaluation of different approaches to annealing schedules and other method optimisations are given by Ingber (1993). As a tool for operational research SA is presented by Eglese (1990), while Koulamas et al. (1994) provide a complete survey of SA applications to operations research problems.

The use of the Simulated Annealing algorithm appears to be a promising approach. We believe that this algorithm could be applied in conjunction with a fast decomposition algorithm to solve efficiently and accurately the aforementioned optimisation problems in much longer production lines.

The remainder of the paper is organised as follows: we first describe the production line model and the problem of our interest followed by the methodology of our approach namely: the performance model, the expansion method used for evaluating the line performance, an overview of the combinatorial optimisation methods, the simulated annealing optimisation method, and the complete enumeration method; we then describe our experimental methodology and present an overview of the numerical and performance results for short and long production lines. In the Appendix following the concluding section, we provide a full tabulated set of the experimental results.

The production line model and the optimal design problem

We define an asynchronous line as one in which every workstation can pass parts on when its processing is complete as long as buffer space is available. Such a line is subject to manufacturing starving and blocking. We assume that the first station is never starved and the last station is never blocked. Therefore we can say that the line operates in a push mode: parts are always available when needed at the first workstation and space is always available at the last workstation to dispose of completed parts.

[line] An \(N\)-workstation multi-machine production line with \(N-1\) intermediate buffers

An \(N\)-station line consists of \(N\) workstations in series, labelled \(M_{1}, M_{2}, \ldots, M_{N}\) and \(N-1\) locations for buffers, labelled \(B_{2}, B_{3}, \ldots, B_{N}\), is illustrated in figure [line]. Each station \(i\) has \(s_{i}\) servers operating in parallel. The buffer capacities of the intermediate buffers \(B_{i}\), \(i=2,\ldots,N\), are denoted by \(q_{i}\), whereas the mean service times of the \(i\) stations (\(i=1,\ldots,N\)) are denoted by \(w_{i}\).

The main performance measure of the production line is the mean throughput, denoted by \(R(\underline{q},\underline{s},\underline{w})\), where \(\underline{q}\)\(=(q_{2},q_{3},\ldots,q_{N})\), \(\underline{s}\)\(=(s_{1},s_{2},\ldots,s_{N})\) and \(\underline{w}\)\(=(w_{1},w_{2},\ldots,w_{N})\).

If \(Q\) denotes the total number of available buffer slots to be allocated to the \(N-1\) buffers and \(S\) the total number of available servers (machines) to allocated to the \(N\) stations then the general version of the optimisation model (first reported by Hillier and So, 1995) is:

\[\max R(\underline{q},\underline{s},\underline{w}) \label{1}\]

subject to:

\[\begin{aligned} \sum_{i=2}^{N} q_{i} = Q, \label{2} \\ \sum_{i=1}^{N} s_{i} = S, \label{3} \\ \sum_{i=1}^{N} w_{i} = N, \label{4}\end{aligned}\]

\[q_{i} \geq 0 \wedge \mathrm{integer} (i=2,3,\ldots,N),\] \[s_{i} \geq 1 \wedge \mathrm{integer} (i=1,2,\ldots,N),\] \[w_{i} > 0 (i=1,2,\ldots,N),\] where \(Q,S\) and \(N\) are fixed constants and \(\underline{q}\), \(\underline{s}\) and \(\underline{w}\) are the decision vectors. The third constraint indicates that the sum of the expected service times is a fixed constant, which by normalisation can be equal to \(N\).

The objective function of throughput, \(R(\underline{q},\underline{s},\underline{w})\), is not the only performance measure of interest. The average WIP, the flow time, the cycle time, the system utilisation, the average queue lengths and other measures are equally important performance measures. However, throughput is the most commonly used performance measure in the international literature.

Optimal allocation methodology

Performance models

The queuing model \(M/M/C/K\) that we use refers to a queuing system where:

• it is assumed that the arrivals are distributed according to the Poisson distribution (or equivalently that the intermediate times between two successive arrivals are exponentially distributed),

• the service (processing) times follow the exponential distribution,

• there are \(C\) parallel servers (machines which are identical at each workstation), and

• the total capacity of the system is finite and equal to \(K\).

While our focus in this paper is on \(M/M/C/K\) approximations for open queuing networks of series-parallel topologies, we also briefly discuss some of the available approaches used for modelling \(M/M/1/K\) systems since most of the literature has focused on \(M/M/1/K\) systems.

Both open and closed systems have been studied by exact analysis although results have been limited. Exact analyses of open two, three, and four node-server models with exponential service are limited by the explosive growth of the Markov chain models for analysing these systems. The analysis of very large Markov chain models has led to effective aggregation techniques for these models Schweitzer and Altiok (1989), Takahashi (1989) but the computation time and power required for these exact results leaves open the need for approximation techniques.

Van Dijk and his co-authors (1988; 1989) have developed some bounding methodologies for both \(M/M/1/K\) and \(M/M/C/K\) systems and have demonstrated their usefulness in the design of small queuing networks. Of course, when doing optimisation of medium and long queuing networks, bounds can be far off the optimum, so robust approximation techniques close to the optimal performance measures are most desirable.

Most approximation techniques appearing in the literature rely on decomposition/aggregation methods to approximate performance measures. One and two node decompositions of the network have been carried out, all with varying degree of success.

The few approximation approaches available in the literature can be classified as follows: Isolation methods, Repeated Trials, Node-by-node decomposition, and Expansion methods. In the Isolation method, the network is subdivided into smaller subnetworks and then studied in isolation Labetoulle and Pujolle (1980), Boxma and Konheim (1981). This method was used by Kuehn (1979) and Gelenbe and Pujolle (1976) but they failed to consider networks with finite capacity.

Closely related to the Isolation method is the Repeated Trials Method, a class of techniques based upon repeatedly attempting to send blocked customers to a queue causing the blocking Caseau and Pujolle (1979), Fredericks and Reisner (1979), Fredericks (1980).

In Node-by-node decomposition, the network is broken down into single, pairs, and triplets of nodes with augmented service and arrival parameters which are then studied separately Hillier and Boling (1967), Takahashi, Miyahara, and Hasegawa (1980), Altiok (1982), Altiok and Perros (1986), Perros and Altiok (1986), Brandwajn and Jow (1988). More general service time approximations appear in the work by Gun and Makowski (1989). The Expansion method is the approach argued for in this paper for computing the performance measures of \(M/M/C/K\) finite queuing networks Kerbache (1984), Kerbache and Smith (1987), Kerbache and Smith (1988). It can be characterised conceptually as a combination of Repeated Trials and Node-by-Node Decomposition where the key difference is that a ‘holding’ node is added to the network to register blocked customers. The addition of the holding node ‘expands’ the network. This approach transforms the queuing network into an equivalent Jackson network which is then decomposed allowing for each node to be solved independently. We have successfully used the Expansion Method to model \(M/M/1/K\) Kerbache and Smith (1988), \(M/M/C/K\) Jain and Smith (1994), Han and Smith (1992), \(M/G/1/K\) Kerbache and Smith (1987), and most recently \(M/G/C/C\) Cheah and Smith (1994), Smith (1991), Smith (1994) queues and queuing networks. In addition, we have also used our Expansion Methodology to model routing Daskalaki and Smith (1989), Daskalaki and Smith (1986), Gosavi and Smith (1990) and optimal resource allocation problems Smith and Daskalaki (1988), Smith (1991), Smith and Chikhale (1995), Singh and Smith (1997).

Expansion method

The Expansion Method is a robust and effective approximation technique developed by Kerbache and Smith (1987). As described in previous papers, this method is characterised as a combination of Repeated Trials and Node-by-node Decomposition solution procedures. Methodologies for computing performance measures for a finite queuing network use primarily the following two kinds of blocking:

Type I:

The upstream node \(i\) gets blocked if the service on a customer is completed but it cannot move downstream due to the queue at the downstream node \(j\) being full. This is sometimes referred to as Blocking After Service (BAS) Onvural (1990).

Type II:

The upstream node is blocked when the downstream node becomes saturated and service must be suspended on the upstream customer regardless of whether service is completed or not. This is sometimes referred to as Blocking Before Service (BBS) Onvural (1990).

The Expansion Method uses Type I blocking, which is prevalent in most production and manufacturing, transportation and other similar systems.

Consider a single node with finite capacity \(K\) (including service). This node essentially oscillates between two states — the saturated phase and the unsaturated phase. In the unsaturated phase, node \(j\) has at most \( K- 1 \) customers (in service or in the queue). On the other hand, when the node is saturated no more customers can join the queue. Refer to figure [Blocking] for a graphical representation of the two scenarios.

Type I Blocking in Finite Queues

The Expansion Method consists of the following three stages :

• Stage I: Network Reconfiguration.

• Stage II: Parameter Estimation.

• Stage III: Feedback Elimination.

The following notation defined by Kerbache and Smith (1987; 1988) shall be used in further discussion regarding this methodology :

\( h:=\) The holding node established in the Expansion method.

\(\Lambda \) := External Poisson arrival rate to the network.

\( \lambda_{j} \) := Poisson arrival rate to node \(j\).

\(\tilde{\lambda_{j}}\):= Effective arrival rate to node \(j\).

\( \mu_{j}\) := Exponential mean service rate at node \( j \).

\(\tilde{\mu_{j}}\) := Effective service rate at node \( j \) due to blocking.

\( p_{K}\):= Blocking probability of finite queue of size \(K\).

\( p_{K}'\):=Feedback blocking probability in the Expansion method.

\( p_{0}^{j}\) :=Unconditional probability that there is no customer in the service channel at node \( j \) (either being served or being held after service).

\(R \) := Mean throughput rate.

\( S_{j}\) := Service capacity (buffer) at node \(j\), i.e. \( S= K \) for a single queue.

Stage I: network reconfiguration

Using the concept of two phases at node \(j\), an artificial node \(h\) is added for each finite node in the network to register blocked customers. Figure [Blocking] shows the additional delay, caused to customers trying to join the queue at node \(j\) when it is full, with probability \(p_{K}\). The customers successfully join queue \(j\) with a probability \((1 - p_{K})\). Introduction of an artificial node also dictates the addition of new arcs with \(p_{K}\) and \((1 - p_{K})\) as the routing probabilities.

The blocked customer proceeds to the finite queue with probability \((1- p_{K}')\) once again after incurring a delay at the artificial node. If the queue is still full, it is re-routed with probability \( p_{K}'\) to the artificial node where it incurs another delay. This process continues till it finds a space in the finite queue. A feedback arc is used to model the repeated delays. The artificial node is modelled as an \(M/M/\infty\) queue. The infinite number of servers is used simply to serve the blocked customer a delay time without queuing.

Stage II: parameter estimation

This stage essentially estimates the parameters \(p_{K}\), \(p_{K}'\) and \(\mu_{h}\) utilising known results for the \(M/M/C/K\) model.

• \(p_{K}\) : Analytical results from the \(M/M/C/K\) model provide the following expression for \(p_{K}\):

\[p_{K} = {{1}\over{c^{K-c} c!}} {\left(\frac{\lambda}{\mu}\right)}^{K} p_{0}\]

where for \((\lambda/c\mu \not = 1)\)

\[p_{0}= \left[ \sum^{c-1}_{n=0} \frac{1}{n!} \left( {\frac{\lambda}{\mu}}\right)^{n} + \frac{(\lambda/\mu)^{c}}{c!} \frac{1 - (\lambda/c\mu)^{K-c+1}}{1 - \lambda/c\mu} \right]^{-1}\]

and for \((\lambda/c\mu = 1)\),

\[p_{0} = \left[\sum_{n=0}^{c-1} \frac{1}{n!} \left(\frac{\lambda}{\mu}\right)^{n} + \frac{(\lambda/\mu)^{c}}{c!}(K- c + 1) \right]^{-1}\]

• \(p_{K}'\): Since there is no closed form solution for this quantity an approximation is used given by Labetoulle and Pujolle obtained using diffusion techniques Labetoulle and Pujolle (1980).

\[p_{K}' = \left[\frac{\mu_{j} + \mu_{h}}{\mu_{h}} - \frac{\lambda [(r_{2}^{K} - r_{1}^{K}) - (r_{2}^{K-1} - r_{1}^{K-1})]} {\mu_{h}[(r_{2}^{K+1} - r_{1}^{K+1}) - (r_{2}^{K} - r_{1}^{K})]}\right]^{-1}\]

where \(r_{1}\) and \(r_{2}\) are the roots to the polynomial:

\[\lambda - (\lambda + \mu_{h} + \mu_{j})x + \mu_{h}x^{2} = 0\]

while, \(\lambda = \lambda_{j} - \lambda_{h}(1 - p_{K}')\) and \(\lambda_{j}\) and \(\lambda_{h}\) are the actual arrival rates to the finite and artificial holding nodes respectively.

In fact, \(\lambda_{j}\) the arrival rate to the finite node is given by: \[\lambda_{j} = \tilde{\lambda_{i}}(1 - p_{K}) = \tilde{\lambda_{i}} - \lambda_{h}\]

Let us examine the following argument to determine the service time at the artificial node. If an arriving customer is blocked, the queue is full and thus a customer is being serviced, so the arriving customer to the holding node has to remain in service at the artificial holding node for the remaining service time interval of the customer in service. The delay distribution of a blocked customer at the holding node has the same distribution as the remaining service time of the customer being serviced at the node doing the blocking. Using renewal theory, one can show that the remaining service time distribution has the following rate \(\mu_h\):

\[\mu_{h} = \frac{2\mu_{j}}{1 + \sigma^{2}\mu_{j}^{2}}\]

where, \(\sigma^{2}\) is the service time variance given by Kleinrock (1975). Notice that if the service time distribution at the finite queue doing the blocking is exponential with rate \(\mu_j\), then: \[\mu_{h} = \mu_j\] the service time at the artificial node is also exponentially distributed with rate \(\mu_j\).

Stage III: feedback elimination

Due to the feedback loop around the holding node, there are strong dependencies in the arrival processes. Elimination of these dependencies requires reconfiguration of the holding node which is accomplished by recomputing the service time at the node and removing the feedback arc. The new service rate is given by:

\[\mu_{h}' = (1 - p_{K}')\mu_{h}\]

The probabilities of being in any of the two phases (saturated or unsaturated) are \(p_{K}\) and \( (1 - p_{K}) \). The mean service time at a node i preceding the finite node is \(\mu_{i}^{-1}\) when in the unsaturated phase and \((\mu_{i}^{-1} + \mu_{h}'^{-1}) \) in the saturated phase. Thus, on an average, the mean service time at the node i preceding a finite node is given by:

\[\tilde{\mu_{i}^{-1}} = \mu_{i}^{-1} + p_{K}\mu_{h}'^{-1}\]

Similar equations can be established with respect to each of the finite nodes. Ultimately, we have simultaneous non-linear equations in variables \(p_{K}\), \(p_{K}'\), \(\mu_{h}^{-1}\) along with auxiliary variables such as \(\mu_{j}\) and \( \tilde{\lambda_{i}} \). Solving these equations simultaneously we can compute all the performance measures of the network.

\[\begin{aligned} \lambda & = & \lambda_{j} - \lambda_{h}(1 - p_{K}') \label{exp:first} \\ \lambda_{j} & = & \tilde{\lambda_{i}}(1 - p_{K}) \label{exp:second} \\ \lambda_{j} & = & \tilde{\lambda_{i}} - \lambda_{h} \label{exp:third} \\ p_{K}' & = & \left[\frac{\mu_{j} + \mu_{h}}{\mu_{h}} - \frac{\lambda[(r_{2}^{K} - r_{1}^{K}) - (r_{2}^{K-1} - r_{1}^{K-1})]} {\mu_{h}[(r_{2}^{K+1} - r_{1}^{K+1}) - (r_{2}^{K} - r_{1}^{K})]}\right]^{-1} \label{exp:fourth} \\ z & = & (\lambda + 2\mu_{h})^{2} - 4\lambda\mu_{h} \label{exp:fifth}\\ r_{1} & = & \frac{[(\lambda + 2 \mu_{h}) - z^{\frac{1}{2}}]}{2\mu_{h}} \label{exp:sixth} \\ r_{2} & = & \frac{[(\lambda + 2 \mu_{h}) + z^{\frac{1}{2}}]}{2\mu_{h}} \label{exp:seventh} \\ p_{K} & = & \frac{1}{c^{K-c} c!} {\left(\frac{\lambda}{\mu}\right)}^{K} p_{0} \label{exp:eigth} \end{aligned}\]

Equations [exp:first] to [exp:fourth] are related to the arrivals and feedback in the holding node. The equations [exp:fifth] to [exp:seventh] are used for solving equation [exp:fourth] with \(z\) used as a dummy parameter for simplicity of the solution. Lastly, equation [exp:eigth] gives the approximation to the blocking probability derived from the exact model for the \(M/M/C/K\) queue. Hence, we essentially have five equations to solve, viz. [exp:first] to [exp:fourth] and [exp:eigth]. To recapitulate, we first expand the network with an artificial holding node; this stage is then followed by the approximation of the routing probabilities, due to blocking, and the service delay in the holding node; and, finally, the feedback arc at the holding node is eliminated. Once these three stages are complete, we have an expanded network which can then be used to compute the performance measures for the original network. As a decomposition technique this approach allows successive addition of a holding node for every finite node, estimation of the parameters and subsequent elimination of the holding node.

Combinatorial optimisation models

The Buffer Allocation Problem (BAP) is perhaps best formulated as a non-linear multiple-objective programming problem where the decision variables are the integers. Not only is the BAP a difficult \(\cal NP-\)hard combinatorial optimisation problem, it is made all the more difficult by the fact that the objective function is not obtainable in closed form to interrelate the integer decision variables \(\bar x\) and the performance measures such as throughput \(R\), work-in-process \(\bar L\), total buffers allocated \(\sum_ix_i\), and other system performance measures such as system utilisation \(\bar \rho\) for any but the most trivial situations. Combinatorial Optimisation approaches for solving problems like the BAP are generally classified as either exact optimal approaches or heuristic ones.

Exact approaches are appropriate for solving small problem instances or for problems with special structure e.g. the Travelling Salesman Problem, which admit optimal solutions. Classical approaches for achieving an optimal solution include Branch-and-Bound, Branch-and-Cut, Dynamic Programming, Exhaustive Search, and related implicit and explicit enumeration methods. The difficulty with utilising these exact approaches for the BAP such as Branch-and-Bound is that the subproblems for which one seeks to compute upper and lower bounds on the objective function are stochastic, non-linear programming problems which are as difficult as the original problem so little is gained by these exact problem decomposition methods.

This dilemma implies that heuristic approaches are the only reasonable methodology for large scale problem instances of the BAP problem. Heuristic approaches can be classified as either classical Non-linear Programming search methods or Metaheuristics.

Non-linear Programming (derivative-free) search Himmelblau (1972) methods such as, to name a few, Hooke-Jeeves, Nelder-Mead simplex methods, PARTAN, Powell’s Conjugate Direction metods, Flexible Tolerance, the Complex Method of Box, and other related techniques have met with varied levels of success in the BAP literature and are viable means of dealing with the BAP because of the non-closed form nature of the non-linear objective function. While many researchers feel that the objective function is concave or pseudo-concave in the decision variables, the discrete nature of the decision variables makes the problem discontinuous and so no derivative information is available.

Metaheuristic methods such as Simulated Annealing, Tabu Search, Genetic Algorithms, and related techniques have not historically been utilised to solve the BAP; in this paper we shall explore the use of Simulated Annealing.

Simulated annealing

Simulated annealing is an optimisation method suitable for combinatorial optimisation 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 (\(\Delta R_{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.

[boltzmann] Probability distribution \(w \sim \exp(\frac{-E}{kT})\) of energy states according to temperature.

Physical matter can be brought into a low-temperature ground state by careful annealing. First the substance is melted, then it is gradually cooled with a lot of time spent at temperatures near the freezing point. If this procedure is not followed the resulting substance may form a glass with not crystalline order and only metastable, consisting of locally optimal structures Kirkpatrick, Gelatt Jr., and Vecchi (1983). During the cooling process the system can escape local minima by moving to a thermal equilibrium of a higher energy potential based on the probabilistic distribution \(w\) of entropy \(S\): \[% Halliday & Resnick (25-15) p. 557 S = k \ln w\] where \(k\) is Boltzmann’s constant and \(w\) the probability that the system will exist in the state it is in relative to all possible states it could be in. Thus given entropy’s relation to energy \(E\) and temperature \(T\): \[% Halliday & Resnick (25-9) p. 552 dS = \frac{d\bar{ }E}{T}\] we arrive at the probabilistic expression \(w\) of energy distribution for a temperature \(T\): \[w \sim \exp(\frac{-E}{kT})\] This so-called Boltzmann probability distribution is illustrated in figure [boltzmann]. The probabilistic ‘uphill’ energy movement that is made possible avoids the entrapment in a local minimum and can provide a globally optimal solution.

The application of the annealing optimisation method to other processes works by repeatedly changing the problem configuration and gradually lowering the temperature until a minimum is reached.

The correspondence between annealing in the physical world and simulated annealing as used for production line optimisation is outlined in table [corr].

Complete enumeration

As a base-rate benchmark for our approach, we used complete enumeration (CE) of all possible buffer and server allocation possibilities. The number of feasible allocations of \(Q\) buffer slots among the \(N-1\) intermediate buffer locations and \(S\) servers among the \(N\) stations increases dramatically with \(Q\), \(S\), and \(N\) and is given by the formula: \[\label{cenum} \begin{array}{l} \left( \begin{array}{c} Q + N - 2 \\ N - 2 \end{array} \right) \left( \begin{array}{c} S + N - 1 \\ N - 1 \end{array} \right) = \\ \frac{(Q+1)(Q+2)\cdots(Q+N-2)} {(N-2)!} \frac{(S+1)(S+2)\cdots(S+N-1)} {(N-1)!} \end{array}\]

All buffer and server combinations can be methodically enumerated by considering a vector \(\underline{p}\) denoting the position within the production line of each buffer or server. Given the vector \(\underline{p}\) we can then easily map \(\underline{p}\) to \(\underline{q}\) or \(\underline{s}\) using the following equation (for the case of \(\underline{q}\)): \[q_{i} = \sum_{j=1}^{Q} \left\{ \begin{array}{ll} 1 & \mathrm{if}\ p_{j} = i \\ 0 & \mathrm{otherwise} \end{array} \right.\] Given then a buffer or server configuration \(\underline{p}\) we advance to the next configuration \(\underline{p}'\) using the function g: \[\underline{p}' = \mathrm{g}(\underline{p}, 1)\] where g is recursively defined (for \(Q\) buffers) as: \[\mathrm{g}(\underline{p}, l) = \left\{ \begin{array}{ll} \begin{array}{l} (p_{1}, \ldots, p_{l}+1, \ldots, p_{Q}) \end{array} & \mathrm{if}\ p_{l} + 1 \leq N \\ \begin{array}{l} (p'_{1}, \ldots, p'_{l-1}, p'_{l+1}, \ldots, p'_{l+1}) \\ \mathrm{where}\ \underline{p}' = \mathrm{g}(\underline{p}, l + 1) \end{array} & \mathrm{otherwise} \end{array} \right.\] Essentially, g maps the vector of positions to a new one representing another line configuration. When the position \(p_{i}\) of buffer or server resource \(i\), as it is incremented, reaches the last place in the line \(N\) (\(p_{i} + 1 = N\)) g is recursively applied to the buffer or server in position \(p_{i+1}\) setting \(p_{i}\)–\(p_{Q}\) to the new value of \(p_{i+1}\). The complete enumeration terminates when all buffers or servers reach the line position \(N\). To enumerate all buffer and server combinations one complete server enumeration is performed for each line buffer configuration.

Experimental design

The main generative procedure we used was Simulated Annealing Kirkpatrick, Gelatt Jr., and Vecchi (1983), Van Laarhoven and Aarts (1987), Spinellis and Papadopoulos (1997) with complete enumeration used where practical in order to validate the SA results. In order to test the approach, 339 cases of small and large production lines from a wide combination of allocation options were searched using SA or CE.

The generalised queuing network throughput evaluation algorithm was initially ported from a vax vms operating system, to a PC-based Intel architecture. Most changes involved the adjustment of numerical constants according to the ieee 488 floating point representation used on the Intel platforms. Subsequently, the algorithm was rewritten in a pure subroutine form so that it could be repeatedly called from within a program and semi-automatically converted from fortran to ansi c.

The SA algorithm used for the production line resource optimisation is given in the figure below.

1. \([\)Set initial line configuration.\(]\) Set \(q_{i} \leftarrow \lfloor Q/N \rfloor\), set \(q_{N/2} \leftarrow q_{N/2} + Q - \sum_{i=2}^{N} \lfloor Q/N \rfloor\), set \(s_{i} \leftarrow \lfloor S/N \rfloor\), set \(s_{N/2} \leftarrow s_{N/2} + S - \sum_{i=1}^{N} \lfloor S/N \rfloor\), set \(w_{i} \leftarrow 1/N\).

2. \([\)Set initial temperature \(T_{0}\).\(]\) Set \(T \leftarrow 0.5\).

3. \([\)Initialize step and success count.\(]\) Set \(S \leftarrow 0\), set \(U \leftarrow 0\).

4. \([\)Create new line with a random redistribution of buffer space, servers, or service rate.\(]\)

   1. \([\)Create a copy of the configuration vectors.\(]\) Set \(\underline{q}' \leftarrow \underline{q}\), set \(\underline{s}' \leftarrow \underline{s}\), set \(\underline{w}' \leftarrow \underline{w}\).

   2. \([\)Determine which vector to modify.\(]\) Set \(R_{v} \leftarrow \lfloor \mathrm{rand} [0 \ldots 2] \rfloor\).

   3. if \(R_{v} = 0\) \([\)Create new line with a random redistribution of buffer space.\(]\) Move \(R_{n}\) space from a source buffer \(q_{R_{s}}\) to a destination buffer \(q_{R_{d}}\): set \(R_{s} \leftarrow \lfloor \mathrm{rand} [2 \ldots N)\rfloor\), set \(R_{d} \leftarrow \lfloor \mathrm{rand} [2 \ldots N)\rfloor\), set \(R_{n} \leftarrow \lfloor \mathrm{rand} [1 \ldots q_{R_{s}}+1)\rfloor\), set \(q_{R_{s}} \leftarrow q_{R_{s}} - R_{n}\), set \(q_{R_{d}} \leftarrow q_{R_{d}} + R_{n}\).

   4. if \(R_{v} = 1\) \([\)Create new line with a random redistribution of server allocation.\(]\) Move \(R_{n}\) servers from source station \(s_{R_{s}}\) to a destination station \(s_{R_{d}}\): set \(R_{s} \leftarrow \lfloor \mathrm{rand} [1 \ldots N)\rfloor\), set \(R_{d} \leftarrow \lfloor \mathrm{rand} [1 \ldots N)\rfloor\), set \(R_{n} \leftarrow \lfloor \mathrm{rand} [1 \ldots s_{R_{s}}+1)\rfloor\), set \(s_{R_{s}} \leftarrow s_{R_{s}} - R_{n}\), set \(s_{R_{d}} \leftarrow s_{R_{d}} + R_{n}\).

   5. if \(R_{v} = 2\) \([\)Create new line with a random redistribution of service rate.\(]\) Move \(R_{n}\) service rate from source station \(s_{R_{s}}\) to a destination station \(s_{R_{d}}\): set \(R_{s} \leftarrow \lfloor \mathrm{rand} [1 \ldots N)\rfloor\), set \(R_{d} \leftarrow \lfloor \mathrm{rand} [1 \ldots N)\rfloor\), set \(R_{n} \leftarrow \mathrm{rand} (0 \ldots w_{R_{s}}]\), set \(w_{R_{s}} \leftarrow w_{R_{s}} - R_{n}\), set \(w_{R_{d}} \leftarrow w_{R_{d}} + R_{n}\).

5. \([\)Calculate energy differential.\(]\) Set \(\Delta E \leftarrow R(\underline{q}, \underline{s}, \underline{w}) - R(\underline{q'}, \underline{s'}, \underline{w'})\).

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 \(\Delta E < 0\) or \(\exp(\frac{-\Delta E}{T}) > \mathrm{rand} (0 \ldots 1)\), set \(\underline{q} \leftarrow \underline{q}'\), set \(\underline{s} \leftarrow \underline{s}'\), set \(\underline{w} \leftarrow \underline{w}'\), set \(U \leftarrow U + 1\).

7. \([\)Repeat for current temperature\(]\). Set \(S \leftarrow S + 1\). If \(S < \) maximum number of steps, go to step 4.

8. \([\)Lower the annealing temperature.\(]\) Set \(T \leftarrow c T\) (\(0 < c < 1\)).

9. \([\)Check if progress has been made.\(]\) If \(U > 0\), go to step 3; otherwise the algorithm terminates. \(\Box\)

The SA procedure was run with the following characteristics based on the number of stations \(N\):

Maximum trials at given temperature

\(100 \times N\)

Maximum successes at given temperature

\(10 \times N\)

Initial temperature

\(0.5\)

Cooling schedule

Exponential: \(T_{K+1} = 0.9 T_{K}\)

Initial line configuration

Equal division of buffers and servers among stations with any remaining resources placed on the station in the middle.

Reported time

Elapsed wall clock time in seconds.

The following facts clarify the use of the evaluation algorithm:

• the throughput used as the line metric and presented in the tabulated results is the throughput at the last station,

• the line topology graph is a linear series of stations allowing for parallel servers, and

• the initial and the effective arrival rate at the first server are set to \(1.5\).

Three batches of tests were planned and executed. One, presented in tables [hr]–[hesb5], was planned in order to compare the results of our approach with those of Hillier and So (1995). A second set, presented in tables [b9]–[be15], was planned in order to compare the approach to the results obtained using a different evaluative procedure Spinellis and Papadopoulos (1997). Finally, the third batch of tests, presented in tables [b]–[bsr], aimed at establishing our method’s efficacy in determining optimal configurations of large production lines. Tables [b9]–[bsr] contain an additional column, \(T\), presenting the program execution time (\(s\)).

Experimental results

The following paragraphs outline the results we obtained by utilising the described methodology on a variety of problems. Where applicable we compare the obtained results with other published data. In our description of the results we sometimes refer to the phenomena of the bowl and the \(L\)-shape.

The bowl phenomenon occurs when for \(j=1,2,\ldots,N\), \[\begin{aligned} w_{j}^{*} & = & w_{N+1-j}^{*} \\ w_{j}^{*} & > & w_{j+1}^{*}~~~\left(1 \leq j \leq \frac{N-1}{2}\right) \\ w_{j}^{*} & < & w_{j+1}^{*}~~~\left(\frac{N-1}{2} \leq j \leq N+1\right)\end{aligned}\] The name arises from the fact that the value of \(w_{1}^{*}=w_{N}^{*}\) is considerably larger than the other \(w_{j}^{*}\), whereas the other \(w_{j}^{*}\) are relatively close, so that plotting the \(w_{j}^{*}\) versus \(j\) gives the shape of the cross-section of the interior of a bowl.

Correspondingly, the \(L\)-phenomenon is the allocation of extra servers to the first station of the production line. This leads to almost the maximum throughput of the line, which is attained when the extra servers are allocated to the last station of the line, having the advantage that it leads to less work-in-process inventory.

Short lines

Tables [hr]–[hesb5] present results for the optimal allocation using SA for lines similar to the ones examined by Hillier and So (1995). In all cases CE results have been calculated for comparison purposes. In comparing the results obtained with those presented by Hillier and So the following important observations can be made:

1. The service rate allocation (table [hr]) does not follow the bowl phenomenon, but diminishes towards the end of the line.

2. The buffer allocation (table [hb]) does not follow the bowl phenomenon, but increases monotonically across the line.

3. The server allocation (table [hs]) for a small number of servers follows a pattern strikingly similar to the one presented in Hillier and So (1995). However, as the number of servers increases, servers tend to accumulate towards the beginning of the line.

4. The server and service rate allocation (table [hsr]) do not exhibit the \(L\)-phenomenon.

5. The buffer and service rate allocation (table [hbr]) results are very similar to the results presented in Hillier and So (1995). As far as the buffer allocation is concerned, buffers tend to accumulate towards the end of the line. The allocation of work however, does not exhibit the bowl phenomenon presented in the other study and follows the usual descending rate across the line.

6. The buffer and server allocation (tables [hsb34], [hsb5]) results are roughly similar to those presented in Hillier and So (1995) in both the server and the buffer allocation vectors. Servers tend to accumulate towards the beginning and middle of the line, whereas buffers tend to accumulate towards the line ends.

7. Finally, the buffer, server, and service rate allocation (table [hbrs]) roughly follows the shape of service rate allocation presented in Hillier and So (1995), but the allocation of buffers and servers is quite dissimilar.

Concerning the behaviour of SA compared to CE the results of the two methods are as follows:

1. Buffer allocation (tables [hb], [heb]) is the same using SA and CE.

2. The server allocation of a large number of servers (tables [hs], [hes]) using SA differs from the allocation using CE but only in the last two experiments, and this is probably due to incorrect results returned by the evaluation algorithm.

3. For buffer and server allocation for up to 4 stations (tables [hsb34], [hesb34]) SA and CE give the same results.

4. For buffer and server allocation of 5 stations (tables [hsb5], [hesb5]) SA and CE differ in the allocated vectors of some configurations, but with only a slight impact on the resulting throughput.

Long lines

The tables ([b15]–[bsr]) present results for optimal allocation using SA for lines of \(N = 4-70\) stations with \(Q=2 N\) buffers and \(S=2 N\) servers. Where practical (i.e. for small lines) CE results have been calculated for comparison purposes. The optimal allocation of a variable number of buffers to a fixed number of servers (9 servers in table [b9]; 15 servers table [b15]) follows a regular pattern: buffers are added to the line from the right to the left. The results of the CE (table [be9]; table [be15]) confirm the SA results. However, this placement strategy is different from the optimal strategy obtained using the decomposition method Dallery and Frein (1993) as the evaluative procedure Spinellis and Papadopoulos (1997).

Performance analysis

The execution time of the SA optimisation algorithm increased exponentially depending on the number of stations (figure [perf]). This increase was a lot better than the combinatorial explosion of CE, and allowed us to produce near-optimal configurations for relatively large production lines in reasonable time. At the extreme end for a 60 station line SA produced a solution for the allocation of 120 buffers, 120 servers, and the service rate configuration in \(4.5\) hours on a 120MHz Pentium-processor system. This solution using SA required the testing of 238,248 different configurations. The dominant factor of the algorithm is the line configuration evaluation taking on our system about 70ms for every configuration test. Multiplying this time factor by the number of all possible buffer and server allocation configurations for the above case — according to equation [cenum] — gives us an approximate evaluation time using CE of \( \left( \begin{array}{c} 120 + 60 - 2 \\ 60 - 2 \end{array} \right) \left( \begin{array}{c} 120 + 60 - 1 \\ 60 - 1 \end{array} \right) 0.07s = 10^{68}s\) or \(10^{62}\) years.

Summary and conclusions

Our approach of using the expansion method for evaluating the production line configurations generated by simulated annealing optimisation allowed us to explore small and large line configurations in bounded execution time. The results obtained with our approach, compared to those of other studies exhibit some interesting similarities but also striking differences between them in the allocation of buffers, numbers of servers, and their service rates. While context dependent, these patterns of allocation are one of the most important insights which emerge in solving very long production lines. The patterns, however, are often counter-intuitive, which underscores the difficulty of the problem we addressed.

Having demonstrated the viability of using an algorithm based on randomisation techniques for optimising large production lines we now plan to explore other optimisation methods based on similar principles such as genetic algorithms to find how they compare to our approach. The simulated annealing algorithm can be fine-tuned in a number of ways. Results for long production lines from different approaches will allow us to tune the algorithm for optimal performance in terms of execution time and derivation of optimal results and propose prescriptive guidelines for implementing production line optimisation systems.

Appendix: tabulated experimental results

[p]

[p]

[p]

[p]

[p]

[p]

[p]

[p]

[p]

[p]