Optimal buffer allocation in finite closed networks with multiple servers

Computers & Operations Research 35 (2008) 2579 – 2598 www.elsevier.com/locate/cor

Optimal buffer allocation in finite closed networks with multiple servers Mustafa Yuzukirmizia,∗ , J. MacGregor Smithb a Department of Industrial Engineering, Kirikkale University, Yahsihan Kirikkale, Turkey b Department of Mechanical and Industrial Engineering, University of Massachusetts Amherst, MA 01003, USA

Available online 22 January 2007

Abstract In this study, we present an optimal buffer allocation procedure for closed queueing networks with finite buffers. The performance measures are evaluated using the expanded mean value analysis, and the solution procedure is incorporated into a nonlinear optimization scheme to arrive at the sub-optimal buffer space vector. The effectiveness of the method is demonstrated through several numerical experiments. Discussions on convergence and computational complexity are also included. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Finite closed queueing networks; Combinatorial optimization

1. Introduction 1.1. Scope and purpose The buffer allocation problem (BAP) is a difficult NP-hard combinatorial optimization problem which involves the integral and effective distribution of buffer space among the workstations in a production line. Most of the research efforts for the optimization of buffer allocation in queueing networks are devoted to open systems. Closed queueing networks have several distinctive features such as the deadlock phenomenon and constant total work-in-process (WIP). In view of that, we present an optimal allocation procedure for finite closed queueing networks which considers the characteristics of these systems. The assumption about the topology of the network is relaxed and the system considered may have series, split and merge configurations and the service stations may also have multiple servers. This appears to be one of the first algorithmic procedures to optimize buffer allocation in closed networks with general topologies and multiple servers. The problem has been formulated as a stochastic, nonlinear, integer programming problem with the objective of optimizing a cost function consisting of the system throughput, average sojourn time and the total buffers allocated. In the solution procedures, we utilize Powell’s unconstrained search algorithm and employ the approximation method of expanded mean value analysis (EMVA) to estimate the performance measures in the closed queueing network. ∗ Corresponding author.

E-mail addresses: [email protected] (M. Yuzukirmizi), [email protected] (J.M. Smith). 0305-0548/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2006.12.008

2580

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

1.2. Design problems Several approaches have been proposed for the systematic solution of the topological network design problems. Smith and Daskalaki [1] decomposed the problem into a set of three interrelated optimization problems based on the scope of the decisions. Although, global optimizations such as by Gavish and Neuman [2] or partially combined optimizations such as by Daskalaki and Smith [3] could be designed, we prefer to extend Smith and Daskalaki’s [1] classification with the addition of the optimal number of customers, which is specific to closed networks. Accordingly, the interrelated optimization problems are: • • • •

optimal topology problem (OTP), optimal routing problem (ORTP), optimal resource allocation problem (ORAP), optimal customer threshold problem (OCTP).

The OTP is concerned with the number and locations of the service stations, identification of their interconnecting edges and general configuration of components. A number of possible mathematical models for the OTPs were studied by Gavish [4], Gavish and Neuman [2], Lee [5] and Ward et al. [6]. The ORTP is concerned with the control of the customer flow through the network. This decision problem arises when there are several different routes that can be used alternatively, and it can be static, quasi-dynamic and dynamic. In the static routing problem, it is assumed that the topology and the network resources are specified. The problem is then concerned with the determination of the optimal routing probabilities. In the quasi-dynamic routing problem, the information about the topology and the traffic rates is periodically updated and the optimal routing probabilities are recomputed in a given interval. In the dynamic routing problem, the time dependent routing vectors are determined and the changes in topology, traffic rates and routing assignments are optimized. Several researchers addressed the problem of the ORTP. Daskalaki and Smith studied the static routing, and then, the dynamic routing problem [7,8], in open finite queueing networks. Tetzlaff [9] has formulated a model for the determination of the optimal route-dependent service rates for closed networks. Tempelmeier and Kuhn [10] presented a systematic review of models for ORTP in their book on flexible manufacturing systems. The ORAP problem is concerned with the optimal allocation of physical resources, i.e. the number of buffers, the number of servers and the related problems of workload distributions. Among several authors, Shanthukumar and Yao [11] first showed that the throughput rate was an increasing concave function of the number of servers, and applied this result to solve the optimal server allocation in closed networks. For a multi-server, two-node, cyclic finite queue, Foley and Park [12] calculated the optimal allocation of fixed number of buffers that maximizes the throughput. The OCTP problem is specific to closed networks and concerned with the threshold number of customers circulating in the system which governs the degree of capacity usage. Kim et al. [13] determined an upper bound on the optimal number of customers with respect to throughput in a serial three-node closed network with finite buffers. Gonzales [14] studied the OCTP problem in single server closed finite networks with general topologies where the problem was formulated as an extension of the BAP. The interrelated problems of optimal routing and optimal buffer allocation in open finite queueing networks were addressed by Daskalaki and Smith [3]. Shanthukumar and Yao [11] studied simultaneous optimal allocation of servers and the optimal number of customers. Foley and Park [12] determined the optimal number of customers, which depends on both the total number of servers and allocation of the waiting spaces in cyclic two-node network. Additionally, various approaches and concepts applicable to optimal network designs can be seen in the volume edited by Smith et al. [15]. This study focuses on the problem of the ORAP, specifically the problem of buffer allocation. It is assumed that the OTOP, the ORTP, and the OCTP problems have already been solved, and the problem of determining how many units of buffer spaces should be assigned to each station has been addressed. 1.3. Buffer allocation problem The BAP is an important design problem for production systems. It is a stochastic, nonlinear and integer programming problem, and essentially concerned with how many buffers should be allocated so that certain objectives are optimized. Allocating large buffers among workstations lessens congestion, decreases the negative effect of blocking and increases

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

2581

the efficiency of a line by smoothing the manufacturing flow. On the other hand, the buffers are costly, may be inadequate, and should be very carefully allocated, otherwise desired performance measures may not be achieved. Therefore, it is essential to study the problem of buffer allocation and design an effective methodology. The BAP is generally constructed as a multi-objective problem and basically has the following mathematical formulation [1,16,17]: Extremize Z = {f1 (x), ¯ ¯ f2 (x), ¯ . . . , fp (x)} s.t. gi (x) ¯ 0, where

 f1 (x) ¯ := Total buffers allocated ( i xi ) f2 (x) ¯ := System Throughput () f3 (x) ¯ := Cycle time ¯ := Work-in-process(WIP) f4 (x) ... fp (x) ¯ := Workstation utilization (i ). gi (x) ¯ := Constraints on the decision variables. In theory, the BAP is an NP-hard combinatorial optimization problem due to its integer programming nature. It is difficult to find a solution to the problem as N and M increase. Also, the objective function is not obtainable in closed form for interrelating the decision variables of buffer sizes with the performance measures of the system, which makes derivative-based approaches problematic. From the viewpoint of the total buffers to be allocated, the BAP can be either constrained or unconstrained. In a constrained problem, the goal is to maximize the achievable production rate with a given total buffer space. In an unconstrained problem, the total buffer space constraint is relaxed and the objective is to maximize a profit function or minimize a cost function. In other perspectives, Gershwin and Schor [18] provided a unified view of the problem types. Researchers select one or more of the objectives suitable to their purpose and attempt to solve their problems at hand with appropriate other objectives as constraints. Several approaches have been suggested as a solution procedure of the BAP, and generally can be categorized as analytical methods, dynamic programming, simulation methods, heuristic search methods and meta-heuristics. A few analytical models have been developed for small networks. In one such study, Foley and Park [12] determined the optimal allocation of a fixed total number of buffers in a two-node cyclic network, where they modeled the two-node network as a generalized semi-Markov process that characterizes the service completion process. In dynamic programming, the problem has been solved in stages with each stage involving exactly one optimizing variable. The computations at the different stages are linked through recursive computations in a manner that yields a feasible optimal solution to the entire problem. Many researchers proposed solution procedures based on dynamic programming. Among others, Huang et al. [19] and Diamantidis and Papadopoulos [20] investigated the buffer allocation strategy of a serially constructed open finite networks utilizing the dynamic programming approach. Simulation methods have the advantage of dealing with a larger class of systems. However, the computation times severely limit their applicability on BAP problems. Hence, they are mostly utilized a few times and prediction models of the performance measures for alternative buffer sizes are implemented. Roser et al. [21] proposed a method based on a single simulation, estimating the effect of additional buffer capacity onto the system. They demonstrated that their method is applicable to large systems, balanced and unbalanced systems and serial and parallel systems. Azadivar and Lee [22] suggested a procedure for a central server closed network using simulation and a stochastic search procedure where the problem was formulated as a constrained stochastic discrete optimization. Heuristic search methods are probably the most common method utilized in finding the optimal or near optimal buffer levels. They tend to resolve the exponential explosion of the number of alternative allocations by sifting quickly through them to discover those which yield close to optimal results. The alternative buffer allocations are iteratively evaluated using either a queueing network model or a stochastic simulation model, in which the system performance measures are obtained. In one of the first research efforts of the BAP formulation, Altiok and Stidham [23] sought to maximize the throughput balanced with the cost from WIP with no upper bound constraint on the total buffers allocated using the search procedure of Hooke and Jeeves [24]. Smith and Daskalaki [1] exploited Powell’s [25] unconstrained search procedure to capture buffer allocation in open finite networks with general topologies. Soon, Singh and Smith

2582

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

[17] extended this approach to the open finite networks to workstations with multiple servers. Smith and Chikdale [26] used Box’s [27] derivative-free sequential search technique in which the initial values are obtained by utilizing Powell’s method. They examined series, merge and split topologies of open finite queueing networks with an overall buffer constraint bound, and investigated how sensitive the buffer constraint was to the optimal buffer allocation. An analysis of other search methods can be found in the book by Ravindran et al. [28]. Metaheuristics are related to search methods and include simulated annealing, tabu search and genetic algorithms. They follow a local improvement paradigm for harnessing the exponential complexity of the solution space. Applications to buffer allocation and a review of related literature are presented using simulated annealing approaches by Spinellis and Papadopoulos [29], genetic algorithms by Bulgak et al. [30] and tabu search by Glover [31], to name a few. Most of the above formulations were employed in open queueing networks. Closed queueing networks have several distinctive features. For instance, total WIP in the system is given and constant, thus, it is trivial to utilize it. Another feature is the deadlock phenomenon, which avoids undersized allocation of buffers. Therefore, any optimization model involving closed finite queueing networks should incorporate these special features. Gonzales [14] successfully developed an optimization method for closed queueing networks which minimizes the cost of total buffers and the cost from the cycle time with a profit from throughput without an upper bound on the total buffers allocated. He also avoided deadlock situations with supplementary tests. The optimization model proposed in this study utilizes the search method of Powell [25]. The EMVA is employed in the evaluation of performance measures such as throughput and the cycle time. However, we will also incorporate and evaluate other methodologies such as enumeration and simulation, where feasible, in assessing the optimization method. 2. Model description We consider a closed finite network, G(M, A) with M service centers with exponential service times and connected with A arcs. The topology of the system is arbitrary and there is a single class of customers with a constant population of N. Accordingly, we shall adopt the following notation: i 1/i Vi rij

:= := := :=

i (j |n) := ci := bi := Ki := PKi :=  :=

Service rate at ith node (i = 1, . . . , M). Mean service time of the jobs at ith node (i = 1, . . . , M). Visit ratio of ith node (i = 1, . . . , M). Routing probability, probability that a job is transferred to the j th node after service completion at the ith node (i  = j, i = 1, . . . , M, j = 1, . . . , M). Conditional probability of having j jobs at the ith node given that there are n jobs in the system (i = 1, . . . , M, j = 0, . . . , n, n = 0, . . . , N ). The number of parallel servers at the ith node (ci 1, i = 1, . . . , M). Waiting space for the ith node (i = 1, . . . , M). Buffer capacity of ith node including the servers Ki = ci + bi (i = 1, . . . , M). Probability that a node i, with finite buffer size of Ki blocks a customer. The cycles in the given closed network.

The performance measures are:  := Overall throughput rate. i := Throughput rate at ith node (i = 1, . . . , M). i := Utilization of ith node (i = 1, . . . , M). Wi :=Average time spent at ith node (i = 1, . . . , M).  W := Cycle time, i.e. average time for a customer to complete one loop or cycle (W = M i=1 Wi Vi ). Qi :=Mean number of customers at ith node (i = 1, . . . , M). Cost and optimization parameters are: C :=The profit from throughput. CW :=The cost of the cycle time.

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

2583

Function Value of Two Node Closed Queue 46 45 44 Function Value

43 42 41 40 39 38 37 36 (3,3)

(3,4)

(3,5)

(4,2)

(4,3)

(4,4)

(4,5)

(5,2)

(5,3)

(5,4)

(5,5)

Buffer Allocation Fig. 1. Objective function value of 2-node finite closed network.

Cx :=The cost of a unit buffer due to the floor space and equipment costs as well as the WIP costs. xi :=The decision variable, i.e. buffer space of node i(xi = bi + ci ). The problem examined in this study is to find the optimal buffer sizes for closed finite queues. Accordingly Minimize s.t.

Z = −C  + CW W + Cx  j ∈

xj N + 1

xi ci

M 

xi

i=1

∀

for i = 1, . . . , M.

The above objective function seeks to maximize the throughput or production rate of the system, and at the same time, minimize the total average waiting delay of a customer from entry to completion. In the formulation of the problem, the buffer sizes, xi , are the decision variables under our optimization control. Since it would be useful to gain some insight into the problem, let us study the basic behavior of the objective function and examine how the decision variables affect it. For this purpose, consider a two-node finite closed network in which the first node has multiple servers with 1 = 0.5, 2 = 1 and c1 = 2, c2 = 1. There are N = 5 customers in the system. A series of simulation experiments was used to obtain the performance measures of the network as the buffer allocation is varied. All simulation runs, consist of 20 replications, each with 11000 time units with 1000 time units of warm-up period, which yielded a small standard deviation in the performance measures. C = 35, CW = 10 and Cx = 1 are the values assigned to the parameters of throughput profit, the cycle-time cost and the buffer cost, respectively. The estimation of the cost parameters is also an important component of the objective function evaluation, as in any optimization method, and discussed in Section 5. Fig. 1 shows the effects of the decision variables on the objective function value. The minimum total number of buffers we can allocate is 6 in order to avoid deadlock. The function value is highly nonlinear as the buffer vector is varied. Furthermore, an optimal minimum is attained with the allocation of K1 =4 and K2 =3, while another suboptimal minimum with K1 = 5 and K2 = 3 exists. In this representation, we note that there is no pre-set value for the total number of buffers to be allocated.

2584

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

3. Expanded mean value analysis (EMVA) The approximation technique that has been proposed in this study uses insights from the expansion method (EM) for modeling series, merging and splitting topologies. The EM’s algorithmic approach has been successful in evaluating performance measures of finite queueing networks. Kerbache and Smith [32] introduced this method to solve a wide range of finite open queueing networks, with exponential service inter-arrival times with varying traffic intensities. EM was further extended to general open queueing networks [33], i.e. general service times and non-exponential i.i.d arrival rates. Soon, this approach was adopted for open networks with M/M/c/K state-dependent service and applied to arbitrarily configured series-parallel network topologies [34]. Moreover, Gonzales [14] successfully adopted the EM to handle closed queueing networks with finite buffers in conjunction with an equalization phase and well-known mean value analysis. The EM is based on a combination of repeated trials and node-by-node decompositions and characterized by the addition of artificial nodes associated with each finite queue to reroute blocked customers to the blocking queues. These additional nodes effectively expand the network and enable the transformation of such systems into equivalent Jackson networks in which each node can be treated independently. In the proposed EMVA approximation, the blocking effect of the closed queue is taken into account through an approximation of the mean effective service time of each server. The mean effective service time is embedded in the state dependent MVA algorithm. Further, marginal queue length probabilities are used to estimate the blocking probabilities along with these effective service rates. Although MVA was derived to provide an exact analysis of product form networks, its appealing algorithmic features, which are numerical stability, capability of providing a reasonable physical interpretation and working directly with the desired statistics, have aroused intense research interest in developing appropriate variants for approximate analysis of non-product form closed queueing networks. 3.1. Queueing network model We consider closed queueing networks with a set of M service centers and A arcs. There is a single class with constant number of N customers. After completing service at center i, a customer, traveling through the network, will next require service at center j with a probability of rij (i  = j, i = 1, . . . , M, j = 1, . . . , M). The service discipline at all-stations are first-come-first-served (FCFS) and service time distributions are exponential. Each station i may have a multiple number, ci (ci = 1, 2, . . .), of servers. Therefore, the station i can serve a maximum number of ci customers, simultaneously (Fig. 2). The service time at station i is 1/i , and the service capacity, the rate at which a station services customers if no blocking occurs is ci i . We assume there are bi waiting spaces for the ith node and customers are blocked after service. This type of blocking arises, when the destination node j is full, at the moment the customer at node i attempts to enter. The customer at node i is blocked and resides in the server, in consequence, preventing the server from beginning service on another customer. We note that, the servers provide spaces for the customers, as well, which leads the number of total buffer spaces to Ki = bi + ci at node i effectively. Servers

Queue

Station Fig. 2. A single station with multiple identical parallel servers.

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

2585

N

j

i Original Network N

i

PK

h

j

1-PK Expanded Network Fig. 3. Expansion of network.

3.2. Effective service rates The method used for estimating the effective service rates ei of node i, is based on the assumptions made in the EM approximation. In the EM, an artificial node h is added for each finite node in the network, effectively expanding the original network (Fig. 3). The name EM alludes to this first stage, which creates a holding node for the overflow of customers due to blocking. After being serviced at node i, a customer attempts to proceed to downstream node j. The customer joins queue j with a probability of 1 − PK , which is the probability of node j having space to accommodate the customer. However, if node j is saturated, the customer is blocked at the previous node i. An additional delay, caused to customers trying to join the queue at j when its full, with probability PK is incurred at the artificial node h. The artificial node is modeled as an M/M/∞ queue. The infinite number of servers is used simply to serve the blocked customer a delay without queueing. In this method, the mean service time of the jobs in the blocked station increases by the remaining service time of the destination node. On average, the mean service time at node i is effectively given by 1 1 1 = + PK . ei i h

(1)

The remaining service time distribution, h , has the same distribution as the remaining service time of the customer being serviced at the node doing blocking and has the following rate: h =

2i+1 1 + 2i+1 2i+1

,

(2)

where 2 is the service time variance given by Kleinrock [37]. If the service time distribution at the finite queue causing the blocking is exponential with j , then h = j . On another note, if the destination node has multiple servers (cj > 1), then the holding service rate would be h = cj j . This comes from the fact that, when a station causes blocking, it is full, and therefore, it is working at the service capacity.

2586

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

0|n

...

1|n

K-1|n

K|n

K|n

Fig. 4. The states for a finite capacity node.

Blocking Probability

K|n

1

...

2n

-K

Possible number of customers being blocked

Fig. 5. Partitioning the blocking probability.

0|n

1|n

...

K-1|n

K|n

K+1|n

...

n|n

Fig. 6. The states for a infinite capacity node in separable networks.

3.3. Blocking probability In this study, we assume production blocking, sometimes referred to as blocking after service (BAS). In the BAS mechanism, a customer upon completion of its service at node i attempts to enter destination node j. If node j at that moment is full, the customer is forced to occupy server i until the destination node becomes available, and node i is blocked. Server i cannot serve any other customer that might be waiting in its queue. In this type of blocking, the blocked server i effectively provides an additional queueing position for customers waiting service at j. The effective queue capacity of blocking node j increases by 1. At this instant, the node j has full capacity of customers, additionally blocks node i. The state change of such a finite node is illustrated in Fig. 4. When there are n (n K) customers in the system, a finite node may have at most K customers. The last state depicts the blocking state where a blocked customer awaits in the upstream node. Moreover, if we map out the state space of a finite network, we will see that the blocking probability of a particular node consists of further partitions. These partitions are the marginal probabilities of having j customers blocked by node i in the upstream nodes. For example, if node j has K customers, it will cause blocking. Since the number of customers in the network is at most n, 1 to (n − K) customers can be blocked, and awaiting for node j to finish. These customers are not necessarily present in the immediate upstream node but can be situated in previous nodes as well. Hence, using these blocking state partitions (Fig. 5), limited space of a finite node can be extended to virtual infinite buffer spaces which represent the possibility of having K + 1, K + 2, . . . , n customers. Concurrently, for an infinite capacity node in a product-form network, the state alterations will be as in Fig. 6. To show the differences from the finite node, the dissimilar states are drawn in dashed lines. With virtual buffer spaces of a finite node, the correspondence to an infinite capacity node can be established. As can be seen in Fig. 6 the partitions of blocking state will be equivalent to K + 1, K + 2, . . . , n customers in the infinite node Assuming this notion, the blocking probability can be expressed as the sum of conditional probabilities of having virtual jobs more than buffer quantity, Ki , at node i when there are n jobs in the equivalent infinite capacity node. This can be formulated as PKi =

n  j =Ki +1

i (j |n) = 1 −

Ki 

i (j |n).

(3)

j =0

Certainly, for a finite capacity node the state changes are not independent from other nodes. However, by introducing effective service rates, the blocking probabilities can be approximated.

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

2587

Table 1 Parameters of the 4-node cyclic network

i ci Ki

I

II

III

IV

4 1 ∞

1 1 5

3 1 ∞

2 1 ∞

Comparison of Blocking Probability in M=4 nodes 1.2 Simulation Approximation BLOKsim BLOKapp BLOKM/M/1/K

1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1

2

3

4

5

6

7

8

9

10

Number of customers Fig. 7. The blocking probability comparison in M = 4 nodes.

We note that while closed queueing networks have been widely studied, to the best of our knowledge, there has been no closed-form expression or other approximation cited for the estimation of the blocking probability for closed finite networks. The blocking probability from the M/M/1/K finite capacity system is often used, and P B i , the probability of finding the finite node i with capacity Ki saturated, is well known as PBi =

i (1 − i )K i

1 − iKi +1

,

(4)

where i is the utilization ratio of station i. We compare our blocking probability estimation calculated by Eq. (3) with expression (4) in a numerical example with the parameters given in Table 1. The network consists of a serially constructed four-node network with single servers. The server at node 1 is the fastest while the one in node 2 is the slowest. We assign a finite space of 5 to node 2 so that it is the only source of blocking and thus, the probability of blocking is high. The other nodes have infinite capacity. The throughput rate from simulation is used for the calculation of utilization ratio, i , which is needed for the M/M/1/K blocking probability. This is from the observation that, the network is cyclic and the throughput rates of all nodes are equal by the conservation of flow. In Fig. 7, the blocking probability estimates for node 2 are plotted on the lower curves of the figure, along with the throughputs in the top two curves from simulation and the approximation algorithm, EMVA, as described in Section 3.4.

2588

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

When there are N = 6 customers in the system, blocking occurs. The probability of blocking rises gradually, as the number of customers increases. Our blocking probability estimation not only provides an exact starting point for the blocking it also mimics the behavior of it. The estimation from M/M/1/K fails to capture the increase in blocking, since the throughput of the system remains steady after N = 6. 3.4. Algorithm The approximation method can be summarized as follows: Step 1: Initialization. For i = 1, . . . , M i (0|0) = 1, Pi (0) = 0. Step 2: Iteration: n = 1, . . . , N: Step 2.1: For i = 1, . . . , M compute the service rates  i (n) =

ni

if nci ,

ci i

if n > ci .

Step 2.2: For i = 1, . . . , M compute the effective service rates ⎧ 1 ⎪ ∀ i for n = 1, ⎪ ⎨ 1 i = ei (n) ⎪ ⎪ ⎩ 1 + Pi+1 (n − 1) 1 ∀ i for n = 2, . . . , N, i (n) h (n) h is the service rate of the holding node, h = i+1 for exponentially distributed service rates for tandem networks. Step 2.3: For i = 1, . . . , M compute the mean response time Wi (n) =

n 

j

e (j ) j =1 i

i (j − 1|n − 1).

Step 2.4: Compute the system throughput n

(n) = M

i=1 Vi Wi (n)

.

Step 2.5: For i = 1, . . . , M compute the conditional probabilities i (j |n) =

(n) i (j − 1|n − 1) ei (j )

i (0|n) = 1 −

n 

for j = 1, . . . , n,

i (j |n).

j =1

Step 2.6: For i = 1, . . . , M compute the blocking probabilities Pi (n) = 1 −

Ki 

i (j |n).

j =0

The algorithm mentioned above has been implemented and was a subject of an extensive study in [35]. Since all the experiments in Section 6 also include a comparison and assessment of the above method, for presentation purposes the exclusive evaluation of the algorithm has been omitted.

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

2589

Fig. 8. Optimization scheme for the buffer allocation.

4. Optimization method Using the analytical technique of EMVA to estimate the throughput and cycle time of the network and employing Powell’s [25] search method, we have established an effective optimization procedure. Similar iterative algorithms have been successfully developed for open queueing networks by Smith and Daskalaki [1], Gonzales [14], Singh and Smith [17]. An illustration of the proposed methodology for the BAP is given in Fig. 8. First, the details of network such as topology, service rates and the number of customers are input. The cost coefficients of the objective function are also entered.

4.1. Initialization After the description of the network and the objective function parameters are input, a proper initial buffer vector is defined. This is particularly advantageous for the nonlinear search method in order to accelerate the convergence of

2590

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

the procedure. Also, this buffer vector is necessary to start the optimization process, and is determined as follows: Step 1: Relax the finite buffer assumption for all nodes in the network, i.e. Ki = N (i = 1, . . . , M), and calculate the performance measures of the network. Step 2: Estimate an appropriate buffer size using the utilization ratios as: 2.1: Calculate the total utilization ratio: T =

M 

i .

i=1

2.2: Calculate the initial buffer size according to: Ki = INT[(1.5 ∗ N ∗ i )/T ]. Step 3: Verify that the initial buffer is deadlock-free. 4.2. Nonlinear unconstrained search The next step is the nonlinear unconstrained search which is the central aspect of the optimization process. We employ Powell’s [25] method which is a classical derivative free method designed to minimize (maximize) a quadratic function f (x) by uni-dimensional searches from an initial point x1 along a set of conjugate directions. The global steps of the algorithm are the following and we follow the exposition by Dixon [36]: Step 4.0: Initialization. 4.1: Define the initial starting buffer vector, x(1) using the above method, or any other initial guess. 4.2: Set  > 0 and choose a set of linearly dependent search directions d(1) , . . . , d(n) . 4.3: Set the iteration number i = 1 and goto Step 5.0. Step 5.0: Buffer vector iterations. 5.1: While i < n define x (i+1) := min f (x (i) + d(i) ). ∗∗ Objective function evaluations, f (.), are done by EMVA. 5.2: Define x (n+3) := 2x (n+1) − x (1) . 5.3: If f (x (n+3) ) > f (x (1) ) then: goto Step 5.0 else x (1) := min f (x (n+1) + (x (n+1) − x (1) )) and establish a new set of directions: d(1) , . . . , d(n) . Step 6.0: Termination. 6.1: If x (n+1) − x (1)  <  then terminate else goto Step 5.2. 6.2: x (1) ← x (n+1) ; i ← 1 and repeat Step 5.0. 4.3. Deadlock tests After the initialization and in the last part of the optimization method, deadlock tests are carried out to ensure the network is deadlock-free. For a deadlock-free network, the total number of buffers in each cycle, i∈ xi , should be one greater than the number of customers, N. In these tests, each cycle is examined. If there is a cycle which violates the theorem, then, the total size of buffers in the cycle is increased by = N + 1 − i∈ xi . The additional buffers are assigned to the stations which gives the best objective value. Steps of the deadlock tests are as follows: Step 1.0: Deadlock-free  condition check. 1.1: For cycle, j , if i∈j xi N + 1 goto Step 3.0, else go to Step 2.0. Step 2.0: Augment evaluations.  2.1: Set = N + 1 − i∈j xi . 2.2: While > 0, Set f ∗ (x) = 10 000 and i ∗ = −1. For ∀i ∈ j , set xi = xi + 1, evaluate f (x) using EMVA. If f (x)f ∗ (x), then set f ∗ (x) = f (x) and i ∗ = i. Reset xi = xi − 1. 2.3: Set xi ∗ = xi ∗ + 1 and = − 1 goto step 2.2.

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

2591

Total Number of Buffers Allocated to Stations 60 55 50 Total Number of Buffers

45 40 35 30 25 20 15 10 5 0 1xE-5 5xE-5 1xE-4 5xE-4 1xE-3 5xE-3 1xE-2 5xE-2 1xE-1 5xE-1

1

(Buffer Cost/Throughput Profit) Ratio Fig. 9. Change of total buffer allocation with the change of Cx /C ratio.

Step 3.0: Termination. 3.1: If all cycles are checked then terminate, else j = j + 1 go to Step 1.0. 5. Cost parameters Powell’s search algorithm is unconstrained. That is, there is no set limit to the total number of buffers to be allocated. However, due to the nature of the modeled system and the optimization procedure, there exist lower and upper threshold limits in the amount of buffers to be assigned. First of all, we have to allocate a sufficient number of buffers to prevent deadlocks, no matter how expensive the buffers are. Moreover, the system reaches its physical limit as we relax the blocking by allocating more buffers and the throughput becomes insensitive to the total number of buffers. These limits also depict how the optimization method will perform with different cost parameters. If buffers are too costly to allocate, compared to the gain from system throughput, then, the total number of buffers allocated by the optimization method is likely close to the deadlock limit. As buffers get cheaper, the total number of allocated buffers will increase. Alternatively, if the buffers are inexpensive, the optimization method allocates sufficient enough number of buffer so that the system reaches its highest achievable throughput. Since any increment in the allocation of more buffers will have an insignificant effect on the objective function, the optimization method will converge and terminate. This is illustrated in Fig. 9. In this example, we evaluate the total number of buffers allocated by the above optimization procedure as we vary the ratio of buffer cost, Cx , to the profit from throughput, C . The system considered is a 3-node network with single servers and equal service rates. The number of customers in the system is N = 15. As the buffers become more expensive we tend to allocate less. Nevertheless, the deadlock tests prevent assigning fewer than total of 16 buffers. Meanwhile, as the ratio (Cx /C ) goes to 0, the total number of buffers increases. Yet, after a certain limit the network reaches its physical capacity and the total number of buffers allocated does not change, no matter how inexpensive the buffers are. 6. Numerical experiments In this section, numerical examples are presented to illustrate the efficacy of the optimization procedure. The experiments are done for the networks with tandem and split-merge topologies with varying service rates, i ,

2592

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

Table 2 Optimal buffer allocation in 2-node network with i = (0.5, 1), ci = (2, 1) N

1 2 3 4 5 8 10 12 15

K

2,1 2,1 2,2 3,3 4,4 7,7 9,9 10,10 14,13

Objective function

Enumeration

Throughput

Cycle time

fEMVA

fs

%

K enum

fe

EMVA

s

%

WEMVA

Ws

%

21.33 15.33 21.00 30.21 40.47 73.78 96.86 120.66 155.24

21.35 15.43 21.14 30.39 40.73 74.04 97.12 119.70 155.14

−0.1 −0.6 −0.7 −0.6 −0.6 −0.4 −0.3 0.8 0.1

2,1 2,1 2,2 3,3 4,3 7,6 9,8 10,10 13,12

21.34 15.42 21.14 30.39 39.79 73.18 96.22 119.70 154.86

0.333 0.6 0.714 0.778 0.818 0.882 0.904 0.906 0.932

0.333(±0.001) 0.598(±0.001) 0.712(±0.002) 0.776(±0.002) 0.815(±0.002) 0.880(±0.003) 0.903(±0.003) 0.911(±0.002) 0.933(±0.003)

−0.1 0.2 0.2 0.3 0.3 0.2 0.2 −0.6 −0.1

3.00 3.33 4.20 5.14 6.11 9.06 11.05 13.23 16.08

3.00(±0.017) 3.33(±0.011) 4.20(±0.012) 5.15(±0.018) 6.12(±0.018) 9.08(±0.039) 11.07(±0.037) 13.16(±0.039) 16.08(±0.067)

0.0 −0.2 −0.2 −0.2 −0.3 −0.2 −0.2 0.6 0.0

Table 3 Optimal buffer allocation in 2-node network with i = (0.2, 0.5), ci = (5, 3) N

5 10 12 15 18 20 22 25

K

5,3 9,6 11,7 14,9 17,10 19,10 21,10 24,10

Objective function

Enumeration

Throughput

Cycle time

fEMVA

fs

%

K enum

fe

EMVA

s

%

WEMVA

Ws

%

53.86 84.19 105.57 139.35 173.41 195.97 218.45 252.01

53.84 84.21 105.82 139.56 173.94 196.12 218.30 251.57

0.1 0.0 −0.2 −0.1 −0.3 −0.1 0.1 0.2

5,3 7,5 9,6 10,9 12,12 13,13 15,13 16,15

53.84 80.86 102.90 136.60 171.13 192.42 214.55 247.72

0.707 0.969 0.983 0.992 0.993 0.991 0.990 0.989

0.708(±0.002) 0.969(±0.004) 0.982(±0.005) 0.991(±0.006) 0.991(±0.005) 0.991(±0.005) 0.991(±0.005) 0.991(±0.005)

0.0 0.0 0.1 0.1 0.2 0.1 −0.1 −0.2

7.06 10.31 12.20 15.10 18.11 20.16 22.21 25.26

7.06(±0.016) 10.31(±0.039) 12.22(±0.065) 15.12(±0.084) 18.16(±0.090) 20.18(±0.100) 22.19(±0.110) 25.22(±0.125)

0.0 0.0 −0.2 −0.1 −0.2 −0.1 0.1 0.2

number of servers, ci , and routing probabilities, pij . For each topology, we have run the optimization method with  We have simulated each different numbers of customers, N, and obtained the suboptimal values of buffer allocation, K. case with the suggested optimal values of the buffer allocation. Arena 5.0 is used as the simulation software and outputs are collected after 20 replications each consisting of 10 000 units run length with 1000 units of warm-up period. The comparison of the cost objective functions, fEMVA and fs , and the throughput of the system EMVA and s are displayed next, as resulted from the optimization method using EMVA and the simulation, respectively. The corresponding relative percentage errors (%) are also provided. Accordingly, the optimal values of K and the resulting objective value, fe , obtained by enumerating the feasible space of buffer allocation are given in relevant tables. To reveal the quality of the optimization method, the cost parameters primarily used in the experiments are C =−35, CW = 10 and Cx = 1. With these parameters, in some instances the optimal buffers are the outcome of the nonlinear unconstrained search, while, the optimization method carried out deadlock tests after the search in the others. We have also included experiments with the parameters of C = −35, CW = 10 and Cx = 0.5 to further examine the method. In assessment of the optimality and the integration of the approximation method of EMVA, the design experiments can be grouped into three sets. 6.1. Experiments with 2-node network We, first, study a small network that consists of two nodes. In Table 2, the stations are balanced in terms of their service rates and the number of servers where i = (0.5, 1) and ci = (2, 1). In Table 3 the service rates of i = (0.2, 0.5) and the number of the servers, ci = (5, 3), are chosen so that the stations are unbalanced.

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

2593

Table 4 Optimal buffer allocation in tandem 3-node network N = 15

N = 20

N = 25

(a) 3-node tandem network with i = (1, 1, 1) and ci = (1, 1, 1) K 10,10,10 172.0 − 171.8(0.1%) fEMVA − fs (%) EMVA − s (%) 0.8697 − 0.8709(−0.1%)  fe ) 10,10,10(172.0) enum.(K,

13,13,13 232.2 − 230.6(0.7%) 0.8914 − 0.8968(−0.6%) 13,13,13(232.2)

16,17,17 292.9 − 290.8(0.7%) 0.9099 − 0.9161(−0.7%) 17,17,16 (292,9)

(b) 3-node tandem network with i = (1, 0.5, 1) and ci = (1, 2, 1) K 9,10,9 171.9 − 171.4(0.3%) fEMVA − fs (%) EMVA − s (%) 0.8619 − 0.864(−0.2%)  fe ) 9,10,9(171.9) enum.(K,

13,13,13 231.7 − 230.5(0.5%) 0.8932 − 0.8973(−0.5%) 13,13,13 (231.7)

16,17,16 292.1 − 290.4(0.6%) 0.9092 − 0.9145(−0.6%) 16,16,16(292,1)

(c) 3-node tandem network with i = (1, 2, 3) and ci = (3, 1, 1) K 7,9,7 32.6 − 30.9(5.4%) fEMVA − fs (%) 1.937 − 1.96(−1.1%) EMVA − s (%)  fe ) 7,9,7(32.6) enum.(K,

9,11,8 63.4 − 60.6(4.6%) 1.938 − 1.969(−1.6%) 9,11,8(63.4)

10,14,8 93.1 − 90.2(3.3%) 1.938 − 1.967(−1.5%) 10,14,8(93.1)

(d) 3-node tandem network with i = (0.333, 0.5, 0.333) and ci = (3, 2, 3) K 8,8,8 171.6 − 169.0(1.5%) fEMVA − fs (%) EMVA − s (%) 0.8464 − 0.8572(−1.3%)  fe ) 8,8,8(171.6) enum.(K,

11,11,11 230.0 − 229.3(0.3%) 0.8783 − 0.8806(−0.3%) 12,11,11(229.9)

16,15,14 289.7 − 289.0(0.2%) 0.9046 − 0.9068(−0.2%) 15,15,15(289.3)

First, we obtained the optimal buffer vector from the optimization algorithm. The networks have been simulated with these parameters and objective function values are compared. We, then, searched for the optimal allocation with enumerating the buffers on simulation. We presented these values as the number of customers in the system, N, is varied. The cost parameters used in the experiments are C = −35, CW = 10 and Cx = 1. Relevant performance measures of the system such as throughput and cycle time are also compared. The 95% confidence half-width of the simulation values through out the numerical examples are given in parenthesis. 6.2. Experiments with 3- and 5-node networks In Tables 4–7, we present numerical experiments with networks that consist of three to five nodes. After obtaining the optimal buffer vector from the algorithm, these networks have been simulated with the obtained buffer capacity values. To test for optimality, exhaustive enumeration has been carried out employing the approximation method of the EMVA. That is, we incremented the feasible buffers systematically, one by one, and performed objective function evaluations using the EMVA. Millions, if not billions, of function evaluations are done and the process took a significant amount of time. We provide a comparison of running times and the number of evaluations in the Section 7. We also provide a comparison of the objective function value and the throughput of the system based on the values obtained from the EMVA and simulation. The cost parameters used in all of the experiments are C = −35, CW = 10 and Cx = 1. The standard deviation of simulations are very small, on the order of [0.0025, 0.0040] for the throughput values for example, and omitted for presentation purposes. We also display the buffer vectors obtained by enumeration and their corresponding objective function values. The first set of experiments are executed with 3-node tandem closed networks. The service rates of stations and the number of servers are varied from identical and balanced to unequal and unbalanced stations (Table 4). Table 5 shows the results with 3-node split closed networks. Note that, merge topologies in these closed networks are virtually the symmetric counterparts. Similar to previous experiments, the service rates and the number of servers are varied from identical to unequal, also considering the routing probabilities. The 5-node experiments has been presented in Table 6 for serially constructed nodes and in Table 7 for the networks with the split topology.

2594

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

Table 5 Optimal buffer allocation in split 3-node network N = 15 (a) 3-node split network with i = (1, 1, 1), ci = (1, 1, 1) and r12 = 0.5

5 K 14  5 141.6 − 142.1(−0.3%) fEMVA − fs (%) EMVA − s (%) 0.9861 − 0.9836(0.3%)  fe ) 14,5,5(141.6) enum.(K, (b) 3-node split network with i = (1, 0.5, 1), ci = (1, 2, 1) and r12 = 0.33

3 K 14  8 143.2 − 144.3(−0.8%) fEMVA − fs (%) EMVA − s (%) 0.983 − 0.9772(0.6%)  fe ) 14,3,8(143.2) enum.(K, (c) 3-node split network with i = (2, 3, 4), ci = (3, 1, 1) and r12 = 0.5

11 K 11  7 −133.3.134.6(−0.9%) fEMVA − fs (%) 5.428 − 5.409(0.4%) EMVA − s (%)  fe ) 11, 11, 7(−133.3) enum.(K,

N = 20

N = 25

5  5 197.4 − 197.9(−0.3%) 0.9857 − 0.9836(0.2%) 19,5,5(197.4)

5  5 253.1 − 253.8(−0.3%) 0.9858 − 0.9836(0.2%) 24,5,5(253.1)

3  9 199.7 − 200.9(−0.6%) 0.9843 − 0.9796(0.5%) 19,3,9(199.7)

4  10 255.7 − 255.4(0.1%) 0.9907 − 0.9915(−0.1%) 24,4,10(255.7)

19

19

15  8 −119.9.119.9(0.0%) 5.568 − 5.567(0.0%) 16, 15, 8(−119.9)

16

24

24

19  9 −104.5.105.4(−0.8%) 5.627 − 5.648(−0.4%) 20, 19, 9(−104.5)

20

Table 6 Optimal buffer allocation in tandem 5-node network N = 20

N = 30

N = 40

= (1, 1, 1, 1, 1) and ci = (1, 1, 1, 1, 1) 7,7,7,7,7 251.5 − 255.9(−1.7%) 0.8161 − 0.8032(1.6%) 7,7,7,7,7(251.5)

10,10,10,10,10 370.5 − 374.7(−1.1%) 0.856 − 0.8466(1.1%) 11,10,10,10,10(370.4)

15,14,14,14,14 490.8 − 492.4(−0.3%) 0.8872 − 0.8843(0.3%) 15,14,14,14,14(490.8)

(b) 5-node tandem network i = (1, 1, 1, 1, 1) with ci = (1, 2, 1, 2, 1) K 9,4,8,3,8 230.9 − 233.7(−1.2%) fEMVA − fs (%) EMVA − s (%) 0.8719 − 0.8626(1.1%)  fe ) 9,4,8,3,8(230.8) enum.(K,

15,5,13,4,13 349.5 − 349.5(0.0%) 0.9058 − 0.9057(0.0%) 15,5,13,4,13(349.5)

22,5,19,5,19 469.4 − 466.8(0.6%) 0.926 − 0.9316(−0.6%) 22,5,19,5,19(469.4)

(c) 5-node tandem network with i K fEMVA − fs (%) EMVA − s (%)  fe ) enum.(K,

= (1, 2, 3, 2, 1) and ci = (3, 2, 1, 2, 3) 7,5,6,4,6 141.8 − 148.6(−4.5%) 1.265 − 1.224(3.4%) 7,5,5,4,6(141.8)

11,7,10,5,10 223.9 − 227.3(−1.5%) 1.320 − 1.304(1.2%) 12,7,10,6,10(223.6)

17,9,15,8,15 309.1 − 306.6(0.8%) 1.365 − 1.375(−0.7%) 17,9,15,7,16(309.1)

(d) 5-node tandem network with i K fEMVA − fs (%) EMVA − s (%)  fe ) enum.(K,

= (1.33, 0.33, 1, 0.67, 0.5) and ci = (1, 5, 2, 2, 3) 6,6,3,4,5 11,9,5,5,9 159.8 − 177.0(−9.8%) 245.2 − 256.9(−4.6%) 1.137 − 1.053(8.0%) 1.207 − 1.160(4.1%) 6,6,3,4,5(159.77) 11,9,5,6,9(245.08)

18,12,7,8,13 334.9 − 339.8(−1.4%) 1.247 − 1.231(1.3%) 18,12,7,8,13(334.854)

(a) 5-node tandem network with i K fEMVA − fs (%) EMVA − s (%)  fe ) enum.(K,

6.3. Experiments with 8-node networks In the experiments for 8-node, closed queueing networks with tandem and split topology have been considered. For these networks enumeration is not feasible to test for optimality. Since it grows exponentially, the amount of alternative buffer allocations is exhaustive. For example, in a network of 8 nodes with single servers and 20 customers, there are 8.917 ∗ 1012 alternative ways to allocate the buffers excluding the ones with deadlock. While the optimality and

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

2595

Table 7 Optimal buffer allocation in split 5-node network N = 30

N = 40

(a) 5-node split network with i = (1, 1, 1, 1, 1), ci = (1, 1, 1, 1, 1) and r12 = 0.33

7 5 K 19  9 8 319.1 − 321.7(−0.8%) fEMVA − fs (%) EMVA − s (%) 0.982 − 0.9747(0.7%)  fe ) 19,7,5,9,8 (319.1) enum.(K,

8 7  10 9 430.3 − 425.5(1.1%) 0.9879 − 0.999(−1.1%) 26,8,7,10,9 (430.3)

(b) 5-node split network with i = (0.5, 1, 1, 1, 1), ci = (3, 2, 1, 1, 2) and r12 = 0.5

5 13 K 13  14 4 200.8 − 208.0(−3.5%) fEMVA − fs (%) EMVA − s (%) 1.4747 − 1.4339(2.8%)  fe ) 14,5,12,13,4 (200.304) enum.(K,

6 16  17 5 280.1 − 291.1(−3.8%) 1.4862 − 1.4371(3.4%) 22,5,14,15,4 (278.93)

(c) 5-node split network with i = (0.5, 1, 1, 0.66, 0.66), ci = (4, 3, 3, 2, 2) and r12 = 0.5

7 6 K 18  12 11 139.1 − 136.0(2.3%) fEMVA − fs (%) 1.9542 − 1.9823(−1.4%) EMVA − s (%)  fe ) 18,7,6,12,11 (139.111) enum.(K,

9 6  12 11 203.9 − 196.8(3.6%) 1.9288 − 1.9798(−2.6%) 24,9,8,13,12 (203.57)

26

19

26

Table 8 Optimal buffer allocation in tandem 8-node network N = 30

N = 40

N = 50

(a) 8-node network with i = (1, 0.5, 1, 0.5, 1, 0.5, 1, 5) and ci = (1, 3, 1, 3, 1, 3, 1, 3) K 6,4,6,4,6,4,6,4 9,5,9,5,9,5,9,5 361.5 − 394.1(−8.3%) 477.2 − 501.0(−4.8%) fEMVA − fs (%) EMVA − s (%) 0.8538 − 0.7864(8.6%) 0.8847 − 0.843(4.9%)

12,7,12,7,12,7,12,6 596.9 − 606.5(−1.6%) 0.9033 − 0.8887(1.6%)

(b) 8-node network with i = (1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6.1, 7) and ci = (1, 1, 1, 1, 1, 1, 1, 1) K 4,12,10,6,4,3,2,1 5,18,15,8,5,3,2,1 303.5 − 334.9(−9.4%) 405.2 − 443.0(−8.5%) fEMVA − fs (%) EMVA − s (%) 0.9501 − 0.871(9.1%) 0.9738 − 0.8971(8.5%)

7,24,18,9,5,3,2,1 509.7 − 550.5(−7.4%) 0.9814 − 0.9126(7.5%)

C = −35, CW = 10, Cx = 0.5 (c) 8-node network with i = (1, 0.6, 0.5, 0.4, 1, 0.6, 0.5, 4) and ci = (1, 2, 3, 4, 1, 2, 3, 4) K 3,6,5,4,3,6,5,4 4,9,7,4,4,9,7,4 317.9 − 353.6(−10.1%) 416.9 − 453.2(−8.0%) fEMVA − fs (%) EMVA − s (%) 0.9047 − 0.8233(9.9%) 0.9395 − 0.8701(8.0%) C = −35, CW = 10, Cx = 0.5

5,13,8,4,5,13,8,4 519.5 − 561.0(−7.4%) 0.9561 − 0.8896(7.5%)

(d) 8-node network with i = (1.72, 1.47, 1.29, 1.45, 1.07, 1.38, 1.27, 1.66) and ci = (1, 1, 1, 1, 1, 1, 1, 1) K 3,2,4,5,8,9,10,4 4,2,4,6,12,13,14,5 279.8 − 301.5(−7.2%) 375.8 − 399.0(−5.8%) fEMVA − fs (%) EMVA − s (%) 1.0234 − 0.9598(6.6%) 1.046 − 0.991(5.5%) C = −35, CW = 10, Cx = 0.5

4,2,5,7,17,16,20,5 474.9 − 494.9(−4.0%) 1.055 − 1.015(3.9%)

sub-optimality remains conjecture on these group, they provide a good perception on the behavior of the buffer allocation. In these sets, we have also included experiments with different cost parameters, where they are explicitly displayed. The results are presented in Table 8 and 9. The presented examples show very good approximations for the underlying methodology. The performance measures generated from the optimization method and the ones obtained from simulation do not differ by more than 10%, in most cases below 4%. The relative percentage errors are consistent and do not increase with the increase of the population,

2596

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

Table 9 Optimal buffer allocation in split-merge 8-node network N = 30

N = 40

(a) 8-node split network with i = (1, 0.5, 0.5, 1, 0.5, 1, 0.5, 1), ci

7 2 3 K 3 16  8 3 1 308.8 − 327.2 − (5.6%) fEMVA − fs (%) EMVA − s (%) 0.9977 − 0.9455(5.5%)

= (1, 3, 3, 1, 3, 1, 3, 1) and r23 = 0.5

8 2 3 5 23  9 3 1 419.5 − 436.9 − (4.0%) 0.9989 − 0.9603(4.0%)

(b) 8-node split network with i = (5, 4, 1, 1, 1, 1.5, 1.5, 1.5), ci = (1, 1, 1, 1, 1, 3, 3, 3) and r23 = 0.25

6 7 8

6 1 1 K 98 21 12  7 4 3  7 14 16 270.6 − 270.3(0.1%) 365.4 − 369.6 − (1.1%) fEMVA − fs (%) 1.1733 − 1.1747(−0.1%) 1.2127 − 1.1992(1.1%) EMVA − s (%) (c) 8-node split network with i = (3.14, 3.18, 1.34, 1.45, 1.86, 1.41, 1.70, 1.92), ci = (1, 1, 1, 1, 1, 1, 1, 1) and r23 = 0.61

5 13 7

5 19 9 15 8 K 11 7  7 5 1  10 7 1 125.6 − 131.9(−4.8%) 188.3 − 194.7(−3.3%) fEMVA − fs (%) EMVA − s (%) 2.0978 − 2.0379(2.9%) 2.1218 − 2.071(2.5%)

Table 10 Computational comparisons in tandem 5-node network Method

Initial K

# of function evaluation

Convergence K

Timea

Above initial. Maximum Arbitrary Enumeration

14,11, 9,14,12 40,40,40,40,40 9,8,7,6,5 –

60 98 93 8.20 ∗ 108

18,12,7,8,13 18,12,7,8,13 18,12,7,8,13 18,12,7,8,13

1s 1s 1s  45 min

a On

Dell Inspiron 600 M with 1.3 GHz.

the system size or the variation of the network topology. We have similar results in the unbiased examples which we randomly generated the service rates of the stations and the routing probabilities (Tables 8(d) and 9(c)). Furthermore, comparison of the enumeration and the optimization procedure indicates that the proposed technique is a very robust search method. In most of the cases, the obtained optimal buffer through the optimization algorithm matches the one from enumeration. Though, there are very slight discrepancies in a few cases, they are very close. One possible reason is that, during the search method integrality of the buffers is relaxed. While passing the values to the function evaluation, they are rounded off. Moreover, the optimization method reveals some counterintuitive buffer allocation strategies. For example, it is better to allocate the buffers in the downstream of the bottleneck station rather than assigning them to this station. This makes sense, because the rate of the bottleneck station is the maximum possible throughput and further blocking of this station will degrade the performance of the system. 7. Computational complexity While the search method is not polynomially bounded in N or M, during the optimization process we have never experienced a convergence problem and the algorithm converged after at most 100 function evaluations including the deadlock tests. In Table 10, we have compared the number of function evaluations and convergence values of the tandem 5-node network from Table 6(d) with N = 40 to demonstrate the efficiency of the method. It is clear that the proposed optimization method provides both the computational efficiency and the solution quality without any compromise. The computational complexity of the algorithm is expected to increase with the addition of more queues, however we conjecture that it would not be extraordinary.

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

2597

8. Conclusion In this article, we have studied the buffer allocation strategy of the closed queueing networks. The assumption about the topology of the network is job-shop-type production systems and the stations may also have multiple servers. An optimal buffer allocation procedure has been proposed which considers the special features of closed queueing network such as deadlock. The solution method involves the minimization of total buffer allocation and average sojourn time while maximizing the system throughput and utilizes the approximation of EMVA. The effectiveness of the method has been demonstrated through several numerical experiments. Although, one limitation of the methodology is that only suboptimal solutions can be guaranteed, we have illustrated that the algorithm finds the optimal value in most cases, and near optimal otherwise. We are confident that using the methodology presented, one can determine a very cost-effective buffer allocation quickly and very precisely. References [1] Smith JM, Daskalaki S. Buffer space allocation in automated assembly lines. Operations Research 1988;36-2:343–58. [2] Gavish B, Neuman I. Capacity and flow assignments in large computer networks. In: Proceedings of IEEE-INFOCOM ’86. 1986. p. 275–84. [3] Daskalaki S, Smith JM. Combining routing and buffer allocation problems in seriesparallel queueing networks. Annals of Operations Research 2004;125:47–68. [4] Gavish B. Formulations and algorithms for the capacitated minimal directed tree problem. Journal of the Association for Computing Machinery 1983;30:118–32. [5] Lee JH. Design and configuration of reconfigurable atm networks with unreliable links. Journal of ETRI 1999;21-4:9–22. [6] Ward J, O’Sullivan M, Shahoumian T, Wilkes Appia J. automatic storage area network fabric design. Conference on file and storage technologies. 2002. p. 203–217. [7] Daskalaki S, Smith JM. The static routing problem in open finite queueing networks. Paper presented at the ORSA/TIMS meeting, Miami, FL, 1986. [8] Daskalaki S, Smith JM. Real time routing in finite queueing networks. In: Perros HG, Altiok T, editors. Queueing networks with blocking. North-Holland: Amsterdam; 1989. p. 313–24. [9] Tetzlaff U. Optimal design of flexible manufacturing systems. In: Tempelmeier H, Kuhn H, editors. Flexible manufacturing systems, vol. 1993. New York: Wiley; 1990. p. 308–14. [10] Tempelmeier H, Kuhn H. Flexible manufacturing systems. New York: Wiley; 1993. [11] Shanthukumar J, Yao D. Optimal server allocation in a system of multi-server stations. Management Science 1987;33-9:1173–80. [12] Foley R, Park B. Optimal allocation of buffers and customers in a two node cyclic network with multiple servers. Operations Research Letters 2002;30:19–24. [13] Kim DS, Kulkarni DM, Lin F. An upper bound for carriers in a three-workstation closed serial production system under production blocking. IEEE Transactions on Automatic Control 2002;47-7:1134–8. [14] Gonzales EA. Optimal resource allocation in closed finite queueing networks with blocking after service, PhD thesis, Department of Mechanical and Industrial Engineering, University of Massachusetts-Amherst, 1997. [15] Smith JM, Gershwin SB, Papadopoulos CT. Performance evaluation and optimization of production lines, vol. 93. Annals of operations research. Baltzer Science Publishers; 2000. [16] Gross D, Harris C. Fundamentals of queueing theory. New York: Wiley; 1985. [17] Singh A, Smith JM. Buffer allocation for an integer nonlinear network design problem. Computers Operations Research 1997;24-5:453–72. [18] Gershwin SB, Schor JE. Efficient algorithms for buffer space allocation. Annals of Operations research 2000;93:117–44. [19] Huang MG, Chang PL, Chou YC. Buffer allocation in flow-shop-type productions systems with general arrival and service patterns. Computer and Operations Research 2000;29:103–21. [20] Diamantidis AC, Papadopoloulos CT. A dynamic programming algorithm for the buffer allocation problem in homogeneous asymptotically reliable production lines. Mathematical problems in Engineering 2004;3:209–23. [21] Roser C, Nakano M, Tanaka M. Buffer allocation model based in a single simulation. In: Chick S, Sanchez Pj, Ferrin D, Morrice DJ, editors. Proceedings of the 2003 winter simulation conference. 2003. p. 1238–46. [22] Azadivar F, Lee Y. Optimum number of buffer spaces in flexible manufacturing systems. In: Papadopoulos HT, Heavey C, Browne J, editors. Queueing theory in man system, vol. 1993, Chapman-Hall; 1986. p. 344–6. [23] Altiok T, Stidham S. The allocation of interstage buffer capacities in production lines. IEEE Transactions 1983;15:251–61. [24] Hooke R, Jeeves TA. Direct search solution of numerical and statistical problems. Journal of the Association for Computing Machinery 1961;8–2:212–29. [25] Himmelblau D. Applied nonlinear programming. New York: McGraw-Hill; 1987 p. 177–180. [26] Smith JM, Chikdale N. Buffer allocation for a class of nonlinear stochastic knapsack problems. Annals of Operations Research 1995;58: 323–60. [27] Box MJ. A new method for constrained optimization and a comparison with other methods. Computer Journal 1965;8:42–52. [28] Ravindran A, Philips DT, Solberg JJ. Operations research—principles and practices. New York: Wiley; 1987. [29] Spinellis DD, Papadopoulos CT. A simulated annealing approach for buffer allocation in reliable production lines. Annals of Operations Research 2000;93:373–84.

2598

M. Yuzukirmizi, J.M. Smith / Computers & Operations Research 35 (2008) 2579 – 2598

[30] Bulgak AA, Diwan PD, Standridge SR. Buffer size optimization in asynchronous assembly systems using genetic algorithms. Computers and Industrial Engineering 1995;28-2:309–22. [31] Glover SB. Tabu search—Part I. ORSA Journal on Computing 1990;1:190–206. [32] Kerbache L, Smith JM. Asymptotic behavior of the expansion method for open finite queueing networks. Computers and Operations Research 1988;15-2:157–69. [33] Kerbache L, Smith JM. The generalized expansion method for open finite queueing networks. European Journal of Operational Research 1987;32:448–61. [34] Jain S, Smith JM. Open finite queueing networks with M/m/c/k parallel servers. Computers and Operations Research 1994;21(3):297–317. [35] Yuzukirmizi M. Finite closed queueing networks with multiple servers and multiple chains. PhD thesis, University of Massachusetts-Amherst, Department of Mechanical and Industrial Engineering, 2005. [36] Dixon LCW. Nonlinear optimisation. London: English University Press; 1972. [37] Kleinrock L, Queueing systems, vol. I: theory. New York: Wiley; 1975.