- Email: [email protected]
Analyzing closed Kanban-controlled assembly systems by iterative aggregation–disaggregation
Computers & Operations Research 26 (1999) 1015}1039
Analyzing closed Kanban-controlled assembly systems by iterative aggregation}disaggregation Jishnu Hazra , Paul J. Schweitzer *, Abraham Seidmann Indian Institute of Management, Bangalore, 560076, India W.E. Simon Graduate School of Business Administration, University of Rochester, Rochester, NY 14627, USA
Abstract In recent years, many companies have begun using WIP-limiting strategies such as CONWIP, Kanban or drum-bu!er-rope to control parts #ows in complex production systems. This paper analyzes assembly systems with a tree structure, random processing times and a constant WIP control system. We present heuristic version of the exact aggregation}disaggregation theory for "nite Markov chains to evaluate the performance of these closed Kanban-controlled assembly systems. Because the approximation is theory based, it provides a framework for further model development, with some possible extensions described in the paper. The approximation has the novel feature of doing simultaneous multiple partitions of the state space, in such a way that the associated aggregate transition rates are mutually consistent. The methodology is a novel approach towards extending aggregation ideas to fork-and-join queueing networks, and it provides several useful operational and analytical insights. Numerical comparisons with simulation show that the proposed approximation computes accurate estimates of the plant throughput. It provides a fast way to assess the performance and economic impact of changes in the total WIP level (or the number of Kanbans), or in the part routes, on the throughput rate of the assembly system. 1999 Elsevier Science Ltd. All rights reserved. Scope and purpose We develop a new methodology for rapid performance analysis of assembly systems with a tree structure, random processing times, and a constant work in process Kanban control system. Extensive evaluations show that this methodology computes accurate estimates of the plant throughput. It provides a fast way for managers to assess the performance and economic impact of changes in the total WIP level (or the number of Kanbans) and in the part routes, on the overall throughput rate of the assembly system. Keywords: Assembly systems; Kanban-controlled; Aggregation; Disaggregation
* Corresponding author. E-mail address: [email protected] (P.J. Schweitzer) 0305-0548/99/$ - see front matter 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 3 0 5 - 0 5 4 8 ( 9 9 ) 0 0 0 2 4 - 6
1016
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1. Introduction Many manufacturing systems require more than one distinct input part. These machines cannot start processing until all input parts from the preceding stages are available at their input bu!er. This results in synchronization delays, which make such queueing systems di$cult to analyze. The multi-stage assembly system is represented by a directed tree as discussed in [1]. We use the generic term part for components and subassemblies. The part in a bu!er of a downstream machine is an assembly of parts from the bu!ers of upstream machines. We model the detailed parts #ows in the tree-structured assembly system as a continuous-time Markov chain (CTMC). We begin, therefore, with the aggregation}disaggregation (A/D) approach to solving CTMCs with large state spaces and then show how to use it for the rapid approximation of complex assembly queueing systems. We expect our results to become an important part of the process analysis activities undertaken by assembly system designers and managers. The need for our research arises from the observation that even complex queueing network models cannot capture the part mating delays, and simulation models of such systems require special e!ort due to the particular part #ows involved. Our results show that the proposed approximation computes accurate estimates of the assembly plant throughput. 1.1. Exact aggregation}disaggregation We consider an ergodic continuous-time Markov chain with "nite state space S and transition rate matrix +r " i, j3S, iOj,. The equilibrium probabilities [n ] (where n "Pr[in state i]), are GH G GZ1 G the solution to the "S" Kolmogoro! equations: n r " nr ; G GH H HG HZ1 G HZ1 G
i3S (discard one),
(1.1a)
n "1. (1.1b) G GZ1 Often "S" is too large to permit direct solution of Eq. (1.1) either by an exact (pivoting) algorithm or even by successive approximations. An alternative numerical approach is Takahashi's aggregation}disaggregation (A/D) procedure [2}4], in which the state space S is partitioned into J groups of `similara states: S"S(1)#S(2)#2#S(a)#2#S(J).
(1.2)
Attention is focused upon the aggregate equilibrium probabilities, n "Pr[in group S(a)]" n , 1)a)J ? G GZ1? and upon the aggregate transition rates,
n G r , 1)a, b)J, aOb R " GH ?@ n IZ1? I HZ1@ GZ1?
(1.3)
(1.4)
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1017
with the interpretation R dt#o(dt)"Pr[in S(b) at time t#dt" in S(a) at time t], aOb. ?@ [R ] inherits the non-negativity and ergodicity properties from [r ]. ?@ GH The A/D approach replaces Eq. (1.1) with the equivalent quadruple set of equations consisting of Eq. (1.4); the aggregation equations: ( ( na R " n R , 1)a)J (discard one) ?@ @ @? @ @ @$?
(1.5a)
@$?
( n "1; ? ? the disaggregation equations for group S(a), 1)a)J: ( n r " n r # n D G, G GH H HG @ @ @ HZ1 G HZ1? G @$? and the disaggregation rate matrix
n H D G" r , @ n HG IZ1@ I HZ1@ with the interpretation
i3S(a)
(1.5b)
(1.6)
(1.7)
D G dt#o(dt)"Pr[in i at time t#dt" in S(b) at time t]. @ Formally, Eq. (1.5) is the set of Kolmogoro! equations for an aggregated J-state CTMC, even though the aggregated process is not Markovian. The quadruple of Eqs. (1.4)}(1.7) is just su$cient to determine the quadruple of unknowns +n , n , R , D G,. The usual iterative A/D algorithm G ? ?@ @ proceeds by alternation between solving Eq. (1.5) for +n , and solving Eq. (1.6) for +n ,. Much ? G smaller sets of linear equations are solved at each step. Upon convergence, Eq. (1.3) provides a numerical check. The A/D framework, however, su!ers from two drawbacks. First, since the size of the problem now is even larger than "S", it is still impractical for problems with enormous state spaces. Second, if the system performance characteristics, such as throughputs, depend only on +n ,, it seems ? cumbersome to have to deal with both +n , and +n , in order to compute just the latter. G ? Ideally one wants a purely aggregated version of the problem, one which involves just the +n , ? and the +R ,. This would overcome both drawbacks. But a glance at the quadruple shows that this ?@ is impossible, because the +n , depend upon the +R ,, which in turn depend on the +n ,. The exact ? ?@ G +n , cannot be computed, therefore, without already knowing +n ,. ? G 1.2. Heuristic aggregation The pursuit of a purely aggregated version, involving just +n , and +R ,, may however, be ? ?@ undertaken as an approximation. One estimates +R , and then solves Eq. (1.5) to obtain an estimate ?@ of +n ,. As the estimate of +R , improves, so does the estimate of +n ,. ? ?@ ?
1018
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
A whole nested family of procedures can be used to estimate +R ,, with increasing accuracy ?@ accompanied by increasing computational e!ort. At one extreme is solving the full quadruple set of equations with perfect accuracy. At the other extreme is a one-pass procedure, in which a single estimate of the conditional probabilities [n / n ] is inserted into Eq. (1.4) to obtain a G IZ1? I single estimate of the +R ,. In between are various heuristic estimates of +R , that may ?@ ?@ involve intermediate models, typically iterative and more coarsely aggregated. These heuristics exploit structural features of the problem at hand, and no general recipe can be given for their construction. 1.3. Multiple partitions This paper describes one such intermediate procedure as applied to closed Kanban-controlled assembly systems. It has the novel feature of doing multiple partitions of the state space, and it generates estimates of +R , by the requirement that the various partitions generate mutually ?@ consistent estimates of the aggregate transition rates. This leads to a xxed-point problem for the +R ,, which itself is solved by iteration. We show that the procedure works, and discuss how it is ?@ capable of even higher accuracy by additional re"nement of the procedure for estimating the +R ,. ?@ The plan of the paper is as follows. Section 2 describes the Kanban system and our notation. Section 3 sets up a single partition of the state space S. Section 4 sets up the multiple partitions and the resulting "xed-point problem. Section 5 describes the computational procedure and results, and Section 6 o!ers further development of the model. Section 7 contains our summary and conclusions.
2. Closed Kanban-controlled assembly system In this paper we study an assembly system operating under CONWIP (Constant Work In Process) control. The goal is to estimate the throughput and the expected local bu!er occupancy for a given topology and a given number of Kanbans. The Kanban system dynamics can be understood from Fig. 1. The assembly system is a directed tree with M machines numbered 1, 2,2, M (here, M"12). Machine f has exponentially distributed assembly times with mean assembly time 1/k , 1)f)M. There is one root node r (here, r"1), a set ¸ of two or more leaf D nodes (here, ¸"+7, 8, 9, 10, 11, 12,), and a set INT of intermediary nodes (here, INT"+2, 3, 4, 5, 6,). New components are loaded at the leaf nodes, processed or assembled at each machine, and forwarded to the machine's output bu!er. For instance, machine 8 takes parts from the input bu!er (whose queue length is n ) and, after completion, forwards them to its output bu!er (whose queue length is n ). That is, n will decrease by one, and n will increase by one. Leaf nodes start the production process of the basic platform (e.g., a computer's motherboard or a car's chassis.) These nodes do not do assembly. Note that machine 5 takes one part from each of its two input bu!ers (whose queue lengths are n and n ) and after assembly forwards them to its output bu!er (whose queue length is n ). Both n and n will decrease by one, and n will increase by one. The root node, machine 1, combines subassemblies coming from nodes 2 and 3. When the root node completes assembly, a fully-assembled part exits the system, and a Kanban (authorization to
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1019
Fig. 1. Part #ows and bu!er structure in a 12-machine assembly system.
produce) as well as raw material is instantaneously and simultaneously delivered to all the leaf nodes. That is, n and n will both decrease by one, while n , n , n , n ,n , and will all increase by one. n In this notation, n denotes the number of parts in the bu!er located downstream from i and GH upstream from j. n includes both parts on queue and parts in the midst of service at j. There is GH su$cient bu!er space that blocking never occurs. Starvation, however, can occur; for example, machine 5 cannot begin assembly unless it has at least one part in each of its input bu!ers, namely, when n *1 and n *1. Aside from the feedback #ow of Kanbans out of the root node, the topology is a directed tree. In general, the tree is unbalanced. Every machine has at least one input bu!er and, aside from the root node, exactly one output bu!er. Technically, this is a closed queueing network with one fork and several joins [5]. The Kanban control ensures that balanced and coordinated #ows of raw materials occur at all the leaf nodes. In open systems, this is often accomplished by MRP systems. A special case of the tree system shown here is the linear assembly line, with one leaf node, which has been analyzed as a CONWIP (constant work-in-process) system. For extensive discussion of models of closed Kanban systems, including CONWIP, see [1,6,7]. Other recent discussions are given in [5,8}11,14]. There is an extensive literature on open assembly systems, but not on closed ones. A convenient modeling viewpoint is that parts never disappear; they are merely mated together. For example a Kanban will simultaneously bring 6 components to the 6 leaf nodes. Four of these six components will subsequently enter the root node from machine 2, while two others will subsequently enter the root node from machine 3. The completed part consists of the 6 assembled components. After exiting the root node, the completed part is `disassembleda into its 6 components, and these 6 components are delivered back to the 6 leaf nodes. The total number of components never changes.
1020
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
The assembly system is started with all bu!ers empty except for the leaf nodes, each of which have their input bu!ers "lled with N components. There will then always be enough components to assemble exactly N "nished goods. We let m denote the number of parts going upstream on a path from machine f to a leaf node. D A convenient consequence of the above assumptions is path independence; m is independent of the D choice of path between f and a leaf node. In Fig. 1, this means: m "n #n "n #n (equality of two paths) m "n #n #n #n "n #n #n #n "n #n #n "etc. (equality of 6 paths)
(2.1)
"N (Note that m "N). P A convenient way of writing the property of path independence is as the recursion m "m #n fOr, where S( f ) denotes the successor to machine f. 1D D D1D The boundary values are m "n if f is a leaf. For Fig. 1, by taking f"8 or 9, D PD m "m #n "m #n . Thus the two paths between machine 5 and the two leaves 8 and 9 have the same number of parts. Path independence is lost if any machine requires more than one part from any of its input bu!ers, or if initial conditions have unequal numbers of parts on the di!erent paths. For simplicity our analysis employs these convenient conditions, but it could be undertaken without them. The network topology is characterized by M"number of nodes, r"root node, INT"set of intermediate nodes, neither root nor leaf, ¸"set of leaf nodes or nodes without predecessors (leaf nodes have only one input bu!er), S( f )"unique successor (downstream) node to node f, fOr [e.g., S(8)"5], P( f )"non-empty set of predecessor (upstream) nodes to node f, f,¸ [e.g., P(1)"+2, 3,]. This topological information, plus N and the +k ,, constitute the full set of parameters for the G problem. The full state space S consists of the set of bu!er occupancy numbers n"+n " 1)i)M, j3S(i),, restricted by the equality of parts on each path and by having N parts GH on each leaf-to-root path. The dynamics can be described as a continuous-time Markov chain (CTMC) with states n"+n , and state space S. The transition rate for machine i (i3INT) to complete an assembly is GH k if n *1 for all j3P(i), G HG 0 otherwise.
If such a transition occurs, the output bu!er level n S(i) increases by one while all input bu!er levels G +n " j3P(i), decrease by one. HG
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1021
The equilibrium probabilities of this CTMC are denoted by n(n), n3S. Since realistic problems can easily have 10 states, it is hopeless to try to calculate them directly from Eq. (1.1). Furthermore, they are not of central interest. The main performance measures of interest are h "throughput of node f, 1)f)M, D calculated from h "k Pr [all input bu!ers for node f are non-empty;] D D U "utilization of node f, 1)f)M, D calculated from
(2.2)
; "h /k ; (2.3) D D D E[n ]"mean occupancy level at the output buwer of node i; 1)i)M, iOr, j"S(i) GH E[n ]"mean occupancy level at the input buwer to each leaf node j, j3L. PH These are calculated from E[n ]" n n(n). GH GH n
(2.4)
The problem can be made dimensionless by measuring the service rates +k , and the throughputs D +h ,, as multiples of k or of min k . D P E E The following properties, Eqs. (2.5) and (2.6), will hold for the exact model and should hold, at least as approximate equalities, for any approximation: h "h "2"h . (2.5) + All machines must have the same throughput, and this cannot exceed min k . Due to WDW+ D starvation, the throughput can be appreciably less than min k . D D The path-independence condition of mean bu!er levels upstream of any machine [e.g., Eq. (2.1)]. (2.6) A more systematic way of writing the path-independence condition is as the recursion Em "Em #En S( f ), all fOr 1D D D plus the boundary values Em "N and Em "En (if f is a leaf). P D PD
(2.7)
3. The partition associated with machine f 3.1. The aggregate states Pick any non-leaf machine f with two or more input bu!ers. For simplicity the description is given as in Fig. 2 for the case in which f has 2 input bu!ers and two predecessors, labeled g and h.
1022
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
Fig. 2. Bu!er content feeding an assembly mode.
[The analysis is simpler if f has one input bu!er, and more cumbersome if f has 3 or more input bu!ers. The framework is the same in all cases.] We generate a partition based on the triple (m , n , n ), of non-negative integers, where D ED FD 0)max[n , n ])m )N. ED FD D That is, a in Eq. (1.2) is this triple a"(m , n , n ), D ED FD and S(a) is the set of states n3S with the given values of m , n , n . All other state information is D ED FD discarded. The aggregate probability [see Eq. (1.3)] that this set of states occurs is denoted by nD(m , n , n ), with the superscript indicating this as a partition based on machine f. S(a) is D ED FD distinguished from S( f ) because a is a triple while f is an integer between 1 and M. It is clear that any accurate analysis of machine f must contain at least information about both of its input bu!er levels n and n . Carrying some additional upstream information via m is useful ED FD D because it supplies m "m !n and hence some partial information about how likely g is to E D ED complete a service, thereby increasing n by one. ED 3.2. Aggregate transition rates According to Eq. (1.5), for the purpose of computing aggregate transition rates, the changes in the triple a"(m , n , n ) behave as a CTMC. The four possible transitions for the triple are as D ED FD follows for f3INT: (a) Machine f has a service completion Next triple b:m , n , n all decrease by one. The aggregate transition rate `R a is D ED FD ?@ k I (n *1 and n *1) and is known. Here I(E) is the indicator function, 1 if E holds and D ED FD 0 otherwise. (b) Machine g has a service completion Next triple b:n increases by one, m and n are unchanged. ED D FD Case 1: g is a leaf. The aggregate transition rate `R a is k I (m !n *1) and is known. ?@ E D ED Case 2: g is not a leaf. The aggregate rate `R a is unknown and will have to be estimated. ?@
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1023
Denote the exact value by kE(m , n , n ). The expression kE(m , n , n ) is to be interpreted D ED FD D ED FD as the e!ective service rate of g as seen by f when in the aggregate state a"(m , n , n ). The D ED FD rate is given by Eq. (1.4) as kE(m , n , n )"k Pr [all of g s input bu!ers are non-empty " m , n , n ] D ED FD E D ED FD In both cases, the rate is zero unless m "m !n *1. E D ED (c) Machine h has a service completion Analogous to (b): if h is not a leaf, there is an unknown e!ective service rate kF(m , n , n ) for D ED FD transition to a new state where n increases by one. FD (d) Kanban arrives simultaneously at all the leaf nodes The next triple b has m increased by one, with n and n unchanged. The aggregate rate D ED FD `R a is unknown and must be estimated. We denote the exact aggregate transition rate `R a ?@ ?@ as cD(m , n , n ). It is to be interpreted as the arrival rate of Kanbans as seen by f while in D ED FD aggregate state a"(m , n , n ). D ED FD The rate is given by Eq. (1.4) as cD(m , n , n )"k Pr[all input bu!ers to root node are non-empty " m , n , n ]. D ED FD P D ED FD Of course this rate vanishes if fOr and m "N, because the root node is starved for parts D coming from machine f. If f3INT, the Kanban arrival rate cD(m , n , n ) is unknown and must be estimated. If f"r, D ED FD the Kanban arrival rate is known exactly: cP(m "N, n , n )"k I(n *1 and n *1). P EP FP P EP FP The aggregate CTMC associated with a non-leaf machine f has as its states the triples (m , n , n ) D ED FD with 0)max[n , n ])m )N, and it involves some unknown aggregate transition rates ED FD D kE, kF and cD. If these rates were known exactly, solution of this aggregate problem (1.5) would give exact values for the aggregate equilibrium probabilities: nD(m , n , n )"Pr[in a state with given values of m , n , n ]. D ED FD D ED FD 3.3. Approximate aggregate transition rates The aggregate transition rates are unknown and must be estimated. We denote estimates with a*, namely kE*, cD*, nD*, h*, Pr[ ]*, En* , and En* . We use approximations of the form (3.1) and D ED FD (3.5) below. Speci"cs about kE* and cD* are given in Section 4.1. The key approximation for kE is kE(m , n , n ) kE*(m !n ); g3INT D ED D ED FD
(3.1)
i.e., the e!ective service rate at g depends only on m . If g is a leaf node, Eq. (3.1) is replaced by the E exact result: kE(m , n , n )"k I(m !n *1), g3¸. D ED FD E D ED
1024
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
We know kE(m , n , n )"0 if m "0, and we therefore require that D ED FD E kE*(m )"0 if m "0. E E
(3.2)
In addition, since kE(m , n , n )"k Pr[g busy " m , n , n ], D ED FD E D ED FD
(3.3)
we know that kE(m , n , n ))k , and we therefore require that D ED FD E kE*(m ))k , m "1, 2,2, N. E E E
(3.4)
The key approximation for cD is cD(m , D cD(m , D
n , ED n , ED
n ) cD*(m ) fOr, D FD n )"k I(n *1 and n *1) f"r, FD P EP FP
(3.5)
i.e., the e!ective Kanban arrival rate seen by fOr depends only on m . D Since cD(n , n , n )"0 when fOr and m "N (the root machine is starved for parts from D ED FD D machine f ), we therefore require that cD*(N)"0, fOr.
(3.6)
cD(m , n , n )"k Pr[r busy " m , n , n ] D ED FD P D ED FD
(3.7)
Since
(The Kanban arrival rate " the completion rate at the root), we know that cD(m , n , n ))k , D ED FD P and we therefore require that cD*(m ))k . D P
(3.8)
3.4. Special cases for f If machine f has only one input bu!er, then the partition is on the pair (m , n ). The analysis is D ED simpler because transitions of type (c) are absent. If f is the root machine r, the analysis is somewhat simpler because m "N, so only n and P EP n enter the CTMC. In this case, transitions (a) and (d) coincide and have rate FP k I(n *1 and n *1) to a new triple where both n and n have decreased by one. P EP FP EP FP 3.5. Model outputs Given these estimated transition rates +kE*, kF*, cD*,, we can solve the aggregate CTMC (1.5) and obtain estimates nD*(m , n , n ) of the exact aggregate equilibrium probabilities D ED FD nD(m , n , n )["`naa]. D ED FD
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1025
3.6. Byproducts of partition at f As by products, we can calculate estimates of throughputs and mean queue lengths from Eqs. (3.9) and (3.10): h*"estimated throughput at f D "k Pr[n *1 and n *1]* [cf. Eq. (2.2)] D ED FD , KD KD "k nD*(m , n , n ) f3¸. (3.9a) D D ED FD KD LED LFD If f is a leaf node, there is no aggregate model such as in Fig. 2 to use for computing h , as in Eq. D (3.9a). However, f is part of the aggregate model
¹"set of predecessors of S( f ) other than f, i.e. siblings of f. associated with S( f ) and we may use h*"h* , f3¸. D 1D En* "estimated queue length at output bu!er of g ED , KD KD " n nD*(m , n , n ). ED D ED FD KD LED LFD
(3.9b)
(3.10)
4. Mutually consistent rates 4.1. Estimating the aggregate rates At this point we have constructed a partition and an aggregate model for one non-leaf machine f, with predecessors g and h, with model inputs kE*, kF* and cD*, and with model outputs nD*(n , n , n ), h*, En* , and En* . Now consider performing a similar partition on every D ED FD D ED FD non-leaf machine. We must provide the model}inputs in a self-consistent way. This will lead to a joint "xed-point problem [Eqs. (4.1)}(4.3) below] for the full set of unknown rates, +kD*(m ), 1)m )N,#+cD*(m ), 0)m )N!1, for all f3INT. D D D D The number of unknowns is 2N"INT""O(NM), and the storage requirements are modest for realistic problems. 4.2. Estimating the ewective service rates The e!ective service rates are estimated as kD*(m )"service completion rate at machine f, given that there are m jobs upstream of f D D "k Pr[n *1 and n *1"m ]*, D ED FD D
1026
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
so kD*(m )"0, D
m "0 D
k KD KD nD*(m , D kD*(m )" D LED LFD D KEDD KFDD nD*(m , L L D
n , ED n , ED
n ) FD , n ) FD
1)m )N D
f3INT, P( f )"+g, h,.
(4.1)
This automatically satis"es Eqs. (3.2) and (3.4).
4.3. Estimating the Kanban arrival rate The Kanban arrival rate, cE*(m ), is estimated by the consistency requirement that g and its E successor f"S(g) `seea the same expected Kanban arrival rate, under the same conditions, namely the same m : E cE*(m )"E[cD*(m ) " m ], 0)m )N, f"S(g) E D E E
0
"
,D EcD (m )Pr[m " m ]* D D E K K *
m "N E fOr, g3INT, 0)m )N!1. E
The last expression is cE*(m )"0, m "N, E E cD*(m )nD*(m , n , n ) D D ED FD , 0)m )N!1, where f"S(g)Or. cE*(m )" E E nD*(m , n , n ) D ED FD
(4.2)
The triple sums are taken over all triples (m , n , n ) which satisfy m )m )N, 0)n )m , D ED FD E D ED D 0)n )m , m !n "m . FD D D ED E When f"r and g3P(r), Eq. (4.2) is replaced by cE*(m )"k Pr[machine r is busy " m ]* E P E "k Pr[n *1 and n *1 " m ]* P EP FP E
0
" k ,\KE nP (m "N, n "N!m , n "k) P I P EP E FP ,\KE nP*(m "N, n "N!m , n "k) I P EP E FP *
m "N, E
(4.3)
0)m )N!1. E
Note that Eq. (4.3) used the exact expression c P(m , n , n )"k I(n *1 and n *1) to evaluate P EP FP P EP FP cE*. Eqs. (4.2) and (4.3) automatically satisfy Eq. (3.6). Eq. (4.3) satis"es (3.8). Induction on the nodes, proceeding upstream from the root, then shows that Eq. (4.2) satis"es (3.8).
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1027
4.4. Overall xxed-point problem Eqs. (4.1)}(4.3) for all f3INT#+r, provide a set of simultaneous relationships among the unknowns +kD*(m ), 1)m )N, f3INT,#+cD*(m ), 0)m )N!1, f3INT,. D D D D They constitute a xxed-point problem for these unknowns. As is the case for other nodal decomposition schemes, this "xed-point set can be solved by successive approximations: kD*(m )"k I(m *1), f3INT#+r,, 0)m )N, D D D D cD*(m )"k I(m (N), f3INT#+r,, 0)m )N. D P D D Loop step: pick next machine f3INT#+r, and generate partition associated with f retrieve the aggregate rates cD* and +kE*, g3P( f ), solve aggregate CTMC for nD* revise kD* [if f3 INT] by Eq. (4.1) revise cE* [if g3P(f) such that g,¸] by Eq. (4.2) or Eq. (4.3) revise h* by Eq. (3.9) D revise En * for all g3P( f ) by Eq. (3.10) ED Termination test: if revised values have not converged, return to loop step. Initialization:
This algorithm typically converges in a few score iterations. The most time-consuming part is the repeated solution of the aggregate CTMCs. We did not explore acceleration by exploiting the skip-free behavior of m , i.e., a generalized birth and death D process. The aggregate CTMCs have O(NR>) states where t"max "P( f )". Even with small t, this D puts a signi"cant restriction on N. When one moves onward from the current choice of f to the next, the nD* can be discarded. Only the +kD*, and +cD*, need be carried forward, so the algorithm's storage requirements are only O(MN). Storing the nD* has the merit of providing a good starting point for the iterative solution of this aggregate CTMC when one returns to it, however, thus reducing computer time. The scheme merits the name iterative aggregation}disaggregation because it alternates between computing the aggregate parameters kD* and cD*, and solving the detailed Kolmogoro! equations for one partition, +nD*(m , n , n ),. Strictly speaking, both of these steps are at the D ED FD aggregate level, but the former is more highly aggregated than the latter. Upon termination, we immediately have the performance measures +h*, and +En* ,. Other D ED performance measures can be computed from nD*.
5. Computational results We start by presenting our results for a 7-machine assembly system with binary tree topology (Fig. 3). It was evaluated both by simulation and by the aggregation}disaggregation algorithm for populations N"6, 9, 12. Table 1 displays the initial numerical results. The principal observations are that the algorithm successfully computes throughput to within 1}2%. It has much lower accuracy (up to 30% error) for predicting expected queue lengths. From a practical system-design perspective, it is important to note that the absolute error in the mean queue length hardly exceeds
1028
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
Fig. 3. The experimental seven-machine assembly system. Every k "1. G
Table 1 Numerical results for test case N
h En "En En "En "En "En En "En "En "En
6
9
12
A/D
Sim
%error
A/D
Sim
%error
A/D
Sim
%error
0.658 2.220 2.600 1.180
0.664 2.360 2.140 1.500
!1.0 !6.0 22.0 !22.0
0.750 3.270 3.900 1.820
0.750 3.490 3.210 2.300
0.0 !3.5 22.0 !20.0
0.804 4.320 5.350 2.430
0.801 4.630 4.270 3.080
0.4 !7.0 23.0 !21.0
Note: h"throughout, En "expected occupancy in output bu!er of machine i. GH %error"100(algorithm!simulation)/simulation. Mean queue lengths are given as average of values at these bu!ers, which should be identical. h"average of the h , which should be identical. D
one part. We attribute this inaccuracy to the inability of the aggregate model to describe correlated behavior among machines. This defect is common to most single-node decomposition methods [12]. This may not be a major practical issue as here we deal with closed queueing systems where bu!er occupancies are also controlled through the size of the population. The assembly system in Fig. 3 has four loops of Kanbans (and parts). We expect, therefore, four relationships of the form En #En #En "N. These held within 1%. The solution violated Eq. (2.5), but the variation among the h* was only about 1%. The average of the h* seems to be D D a good estimate of h. Similarly, the solution violated Eq. (2.7), but the variation among di!erent paths was small as well. Section 6.2 sketches one way to force equalities in (2.5) and (2.7). We "nd that the assembly system throughputs are only 2/3 as large as the ideal value, min +k ,"1, for N"6. That happens because of inherent delays among mating parts as evidD G enced by the increased bu!er occupancies in front of those machines. The results in Table 1 show that the system throughput increases in a concave fashion to 4/5 of the ideal value as we double the number of Kanbans (N"12).
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1029
Table 2 Input data for the seven-machine network in Fig. 1 Data set k
k
k
k
k
k
k
1 2 3 4
1 0.8 0.8 1
1 0.8 1 1
1 0.8 1 0.8
1 0.8 1 0.8
1 0.8 0.8 0.8
1 0.8 0.8 0.8
1 1 1 1
Table 3 Input data for the six-machine network in Fig. 4 Data set k
k
k
k
k
k
1 2 3 4
1 1 1 1
1 0.8 0.8 1
1 1 1 0.8
1 1 1 0.8
1 1 0.8 1
1 1 1 1
Table 4 Input data for the "ve-machine network in Fig. 5 Data set
k
k
k
k
k
1 2 3
1 1 1
1 0.8 1
1 0.8 1
1 1 0.8
1 1 0.8
Table 5 Input data for the "ve-machine network in Fig. 6 Data set
k
k
k
k
k
1 2 3 4
1 1 1 1
1 1 1 1
0.8 0.8 0.8 1
0.8 1 1 1
0.8 0.8 1 1
We conducted extensive computations to better understand the accuracy performance of the algorithm. The results reported here are based on four di!erent network topologies (Figs. 1 and 4}6). For each topology, three or four di!erent data sets are used with three di!erent values of N, the number of Kanbans (N"6, 9 and 12). Kanban #ows from the root node (always numbered 1 in the diagram) to the leaf nodes. The data sets are given by Tables 2}5. Our results in Tables 6}19
1030
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
Table 6 Results for the seven-machine network and Data set 2 in Table 2 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.540 1.977 1.977 2.751 2.751 2.751 2.751 1.275 1.275 1.275 1.275
0.558 2.048 1.968 2.347 2.258 2.338 2.386 1.606 1.694 1.695 1.647
!3.23 !3.47 0.46 17.21 21.83 17.66 15.30 !20.61 !24.73 !24.78 !22.59
0.613 2.803 2.803 4.181 4.181 4.181 4.181 2.008 2.008 2.008 2.008
0.626 2.954 2.735 3.570 3.421 3.576 3.721 2.477 2.625 2.689 2.544
!2.08 !5.11 2.49 17.11 22.22 16.92 12.36 !18.93 !23.50 !25.33 !21.07
0.651 3.602 3.602 5.661 5.661 5.661 5.661 2.727 2.727 2.727 2.727
0.666 3.807 3.501 4.850 4.595 4.860 5.050 3.344 3.598 3.639 3.449
!2.25 !5.38 2.88 16.72 23.20 16.48 12.10 !18.45 !24.21 !25.06 !20.93
Table 7 Results for the seven-machine network and Data set 3 in Table 2 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.587 2.235 1.934 2.804 2.804 2.544 2.544 0.962 0.962 1.523 1.523
0.600 2.373 1.941 2.445 2.389 2.023 2.100 1.182 1.238 2.037 1.959
!2.17 !5.82 !0.36 14.68 17.37 25.75 21.14 !18.61 !22.29 !25.23 !22.26
0.662 3.239 2.760 4.363 4.363 3.753 3.753 1.402 1.402 2.489 2.489
0.670 3.544 2.713 3.808 3.727 2.946 3.103 1.648 1.729 3.341 3.185
!1.19 !8.61 1.73 14.57 17.06 27.39 20.95 !14.93 !18.91 !25.50 !21.85
0.707 4.108 3.628 6.091 6.091 4.849 4.849 1.800 1.800 3.526 3.526
0.708 4.735 3.409 5.218 5.149 3.822 4.079 2.047 2.117 4.768 4.512
!0.14 !13.24 6.42 16.73 18.29 26.87 18.88 !12.07 !14.97 !26.05 !21.85
show that the algorithm computes the system throughput with an error of 5% or less in all cases. That error was not a!ected by the value of N. Most of the estimates for the bu!er occupancy were accurate within 30%, or better and within roughly 1 job in absolute value. Estimation errors for bu!er occupancies were greater for upstream machines than for downstream machines. Again, the magnitude of that estimation error was not necessarily correlated with N. In all cases, increasing N results in a concave increase in the system throughput h. The e!ect of changes in N on h depends of the topology of the assembly system, as can be seen by comparing the two "ve-machine topologies studied. The results of our comparative study are presented in Tables 6}19.
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1031
Table 8 Results for the seven-machine network and Data set 4 in Table 2 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.575 2.063 2.063 2.479 2.479 2.479 2.479 1.460 1.460 1.460 1.460
0.592 2.133 2.074 2.005 1.914 1.960 2.031 1.862 1.953 1.966 1.895
!2.87 !3.28 !0.53 23.64 29.52 26.48 22.06 !21.59 !25.24 !25.74 !22.96
0.650 2.956 2.956 3.656 3.656 3.656 3.656 2.391 2.391 2.391 2.391
0.663 3.072 2.919 2.953 2.741 2.861 3.021 2.975 3.187 3.220 3.061
!1.96 !3.78 1.27 23.81 33.38 27.79 21.02 !19.63 !24.98 !25.75 !21.89
0.694 3.808 3.808 4.799 4.799 4.799 4.799 3.397 3.397 3.397 3.397
0.703 3.925 3.707 3.829 3.552 3.698 3.995 4.246 4.523 4.596 4.298
!1.28 !2.98 2.72 25.33 35.11 29.77 20.13 !20.00 !24.89 !26.09 !20.96
Table 9 Results for the six-machine network of Fig. 4 and Data set 1 in Table 3 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.691 1.910 2.539 2.784 2.784 1.803 1.304 1.304 1.642
0.683 2.216 2.751 2.235 2.175 1.639 1.549 1.609 1.611
1.17 !13.81 !7.71 24.56 28.00 10.01 !15.82 !18.96 1.92
0.775 2.842 3.664 4.161 4.161 2.719 1.996 1.996 2.604
0.766 3.309 4.040 3.357 3.263 2.489 2.335 2.428 2.470
1.17 !14.11 !9.31 23.95 27.52 9.24 !14.52 !17.79 5.43
0.825 3.822 4.699 5.552 5.552 3.636 2.624 2.624 3.656
0.813 4.414 5.371 4.453 4.353 3.319 3.134 3.234 3.311
1.48 !13.41 !12.51 24.68 27.54 9.55 !16.27 !18.86 10.42
6. Further re5nements We sketch several ways that the accuracy of the model can be improved with increased computational e!ort. This demonstrates the #exibility of the A/D framework for developing whole families of approximations. 6.1. More-detailed partitions The current partitioning (m , n , n ) on f and its immediate predecessors, g and h, does not D ED FD provide information either upstream of g and h or downstream of f. A "ner partition would employ
1032
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
Table 10 Results for the six-machine network of Fig. 4 and Data set 2 in Table 3 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.654 2.335 2.006 2.535 2.535 2.498 1.130 1.130 1.483
0.647 2.579 2.154 2.022 1.969 2.401 1.399 1.452 1.445
1.08 !9.46 !6.87 25.37 28.75 4.04 !19.23 !22.18 2.63
0.728 3.695 2.650 3.661 3.661 4.121 1.643 1.643 2.220
0.719 4.131 2.822 2.881 2.798 4.102 1.989 2.072 2.076
1.25 !10.55 !6.09 27.07 30.84 0.46 !17.40 !20.70 6.94
0.768 5.195 3.133 4.718 4.718 5.888 2.087 2.087 2.971
0.755 5.942 3.248 3.579 3.490 6.176 2.479 2.568 2.576
1.72 !12.57 !3.54 31.82 35.19 !4.66 !15.81 !18.73 15.33
Table 11 Results for the six-machine network of Fig. 4 and Data set 3 in Table 3 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.623 2.638 1.587 2.352 2.352 2.270 1.009 1.009 2.133
0.615 2.892 1.743 1.844 1.796 2.125 1.265 1.313 2.133
1.30 !8.78 !8.95 27.55 30.96 6.82 !20.24 !23.15 0.00
0.689 4.292 1.956 3.303 3.303 3.506 1.405 1.405 3.534
0.681 4.721 2.202 2.541 2.455 3.389 1.738 1.824 3.409
1.17 !9.09 !11.17 29.99 34.54 3.45 !19.16 !22.97 3.67
0.722 6.139 2.232 4.136 4.136 4.769 1.725 1.725 4.999
0.714 6.823 2.497 3.060 3.003 4.737 2.117 2.174 4.766
1.12 !10.02 !10.61 35.16 37.73 0.68 !18.52 !20.65 4.89
Table 12 Results for the six-machine network of Fig. 4 and Data set 4 in Table 3 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.625 1.372 3.113 2.817 2.817 1.525 1.811 1.811 1.340
0.627 1.631 3.236 2.232 2.134 1.397 2.138 2.236 1.366
!0.32 !15.88 !3.80 26.21 32.01 9.16 !15.29 !19.01 !1.90
0.691 1.799 4.866 4.158 4.158 2.148 3.042 3.042 1.964
0.696 2.158 5.067 3.347 3.105 1.980 3.495 3.737 1.953
!0.72 !16.64 !3.97 24.23 33.91 8.48 !12.96 !18.60 0.56
0.729 2.179 6.713 5.304 5.304 2.748 4.516 4.516 2.522
0.731 2.548 7.217 4.398 4.045 2.418 5.054 5.408 2.366
!0.27 !14.48 !6.98 20.60 31.12 13.65 !10.65 !16.49 6.59
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1033
Table 13 Results for the "ve-machine network of Fig. 5 and Data set 1 in Table 4 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.733 2.189 2.189 1.968 1.968 1.829 1.829
0.701 2.670 2.534 1.679 1.737 1.651 1.730
4.56 !18.01 !13.61 17.21 13.30 10.78 5.72
0.816 3.153 3.153 2.897 2.897 2.938 2.938
0.781 3.943 3.754 2.569 2.620 2.488 2.627
4.48 !20.04 !16.01 12.77 10.57 18.09 11.84
0.867 4.060 4.060 3.753 3.753 4.179 4.179
0.826 5.230 5.017 3.439 3.461 3.331 3.523
4.96 !22.37 !19.08 9.13 8.44 25.46 18.62
Table 14 Results for the "ve-machine network of Fig. 5 and Data set 2 in Table 4 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.650 2.125 2.125 2.424 2.424 1.437 1.437
0.632 2.448 2.309 2.215 2.297 1.337 1.394
2.85 !13.19 !7.97 9.44 5.53 7.48 3.08
0.715 3.052 3.052 3.866 3.866 2.071 2.071
0.697 3.531 3.325 3.653 3.736 1.816 1.939
2.58 !13.57 !8.21 5.83 3.48 14.04 6.81
0.754 3.923 3.923 5.407 5.407 2.662 2.662
0.731 4.542 4.284 5.234 5.348 2.225 2.368
3.15 !13.63 !8.43 3.31 1.10 19.64 12.42
Table 15 Results for the "ve-machine network of Fig. 5 and Data set 3 in Table 4 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.650 2.004 2.004 1.646 1.646 2.337 2.337
0.631 2.453 2.292 1.358 1.410 2.189 2.298
3.01 !18.30 !12.57 21.21 16.74 6.76 1.70
0.710 2.828 2.828 2.254 2.254 3.906 3.906
0.698 3.572 3.224 1.896 1.978 3.533 3.798
1.72 !20.83 !12.28 18.88 13.95 10.56 2.84
0.744 3.577 3.577 2.701 2.701 5.708 5.708
0.732 4.652 4.117 2.329 2.393 5.019 5.491
1.64 !23.11 !13.12 15.97 12.87 13.73 3.95
1034
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
Table 16 Results for the "ve-machine network of Fig. 6 and Data set 1 in Table 5 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.609 1.574 3.613 2.719 2.719 2.387 1.709 1.709
0.624 1.645 3.652 2.244 2.118 2.348 2.112 2.238
!2.40 !4.32 !1.07 21.17 28.38 1.66 !19.08 !23.64
0.672 2.327 5.047 3.935 3.935 3.952 2.740 2.740
0.688 2.332 5.176 3.301 3.011 3.824 3.367 3.657
!2.33 !0.21 !2.49 19.21 30.69 3.35 !18.62 !25.08
0.706 3.184 6.091 4.945 4.945 5.913 3.870 3.870
0.722 3.012 6.482 4.270 3.793 5.518 4.718 5.196
!2.22 5.71 !6.03 15.81 30.37 7.16 !17.97 !25.52
Table 17 Results for the "ve-machine network of Fig. 6 and Data set 2 in Table 5 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.635 1.756 3.377 3.146 2.361 2.623 1.103 1.881
0.645 1.841 3.447 2.777 1.750 2.553 1.383 2.410
!1.55 !4.62 !2.03 13.29 34.91 2.74 !20.25 !21.95
0.700 2.682 4.556 4.748 3.296 4.444 1.575 3.020
0.707 2.739 4.743 4.363 2.301 4.257 1.898 3.961
!0.99 !2.08 !3.94 8.82 43.24 4.39 !17.02 !23.76
0.733 3.738 5.412 6.112 4.322 6.591 2.142 3.936
0.736 3.656 5.799 6.061 2.666 6.201 2.283 5.678
!0.41 2.24 !6.67 0.84 62.12 6.29 !6.18 !30.68
Table 18 Results for the "ve-machine network of Fig. 6 and Data set 3 in Table 5 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.673 2.090 3.001 2.680 2.680 2.999 1.232 1.232
0.679 2.224 3.107 2.238 2.150 2.893 1.538 1.626
!0.88 !6.03 !3.41 19.75 24.65 3.66 !19.90 !24.23
0.742 3.438 3.697 3.810 3.810 5.303 1.746 1.746
0.744 3.659 3.817 3.168 3.015 5.183 2.173 2.326
!0.27 !6.04 !3.14 20.27 26.37 2.32 !19.65 !24.94
0.782 4.810 4.320 4.962 4.962 7.690 2.225 2.225
0.773 5.470 4.107 3.893 3.677 7.893 2.637 2.853
1.16 !12.07 5.19 27.46 34.95 !2.57 !15.62 !22.01
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1035
Table 19 Results for the "ve-machine network of Fig. 6 and Data Set 4 in Table 5 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.698 1.744 4.094 1.906 2.881 2.881 1.379 1.379
0.710 1.888 4.121 1.879 2.420 2.336 1.691 1.778
! 1.69 !7.63 !0.66 1.44 19.05 23.33 !18.45 !22.44
0.778 2.625 6.044 2.955 4.289 4.289 2.091 2.091
0.789 2.825 6.170 2.830 3.642 3.461 2.533 2.714
!1.39 !7.08 !2.04 4.42 17.76 23.92 !17.45 !22.96
0.820 3.566 7.898 4.101 5.681 5.681 2.750 2.750
0.834 3.788 8.178 3.822 4.810 4.623 3.402 3.589
!1.68 !5.86 !3.42 7.30 18.11 22.89 !19.17 !23.38
Fig. 4. A six-machine network model. (Machines 4}6 are leaf nodes).
additional state variables to capture the additional information. As the partition gets "ner, the model becomes more accurate. 6.2. Adding tunable parameters Estimates (4.1)}(4.3) of the aggregate transition rates can be augmented with additional parameters in order to meet additional conditions on the system. This is illustrated by one approach to ensure that Eqs. (2.5) and (2.7) hold: Expression (4.2) for the estimated Kanban arrival rate cD* as seen by machine f, for f3INT, can be modi"ed by incorporating an additional parameter a . This is justi"ed by the recognition that D
1036
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
Fig. 5. A "ve-machine network model. (Machines 4 and 5 are leaf nodes).
Fig. 6. A six-machine network model. (Machines 3}5 are leaf nodes).
the approximation Eq. (3.5) is not exact. The new parameter helps determine cD* and hence h , and D it should be chosen so that h "h . Satisfaction of this condition for all fOr is Eq. (2.5). D P Similarly, expression (4.1) for the e!ective service rate kD* at machine f3INT can be modi"ed by employing an additional parameter b . This is justi"ed by the recognition that approximation (3.1) D , and it should be chosen so is not exact. This new parameter helps determine kD* and hence En D1D that En "Em !Em . D1D 1D D Satisfaction of this condition for all f3 INT will ensure that the path-independence condition (2.7) holds.
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1037
6.3. Partitions with multiple machines Simulation shows the existence of a `#oating bottlenecka, where completion of an assembly at a slow machine triggers a quick succession of other assembly operations downstream, with this cascade stopping when the platoon of parts reaches a new downstream bottleneck. The aggregate model of f fails to capture this dynamic correlated behavior of several machines, which may be the source of the large errors in computing mean queue lengths, a common failing of `single-nodea decompositions. The appropriate remedy is a partition that simultaneously encompasses several assembly machines. This, unfortunately, leads to an explosion of the number of the states. A coarser partition of the Cartesian-product state space seems appropriate. Computational results for blocking networks [13] suggest that there are diminishing marginal bene"ts from incorporating more than a few machines in the aggregate model. 6.4. Hybrid schemes Hazra [1,7] describes a di!erent aggregation approach for evaluating assembly networks, and surveys other available techniques. Each technique has its own strengths and weaknesses. An open issue is how to construct hybrid approaches that combine the bene"ts of several pre-existing schemes. For example, Hazra's [1] scheme has only 3% error in computing mean queue lengths, while the present A/D scheme has more accurate throughputs. Could the aggregate service rates from one scheme be used as inputs to the other? The goal would be to achieve the advantages of both: only 1% error in throughput and only 3% error in mean queue length. 6.5. Triangular decomposition of the xxed-point problem There is structure in the "xed-point equations that has not yet been exploited. The in#uence diagram for the network in Fig. 1 is given by Fig. 7. Fig. 7 shows the interactions in the "xed-point system, and mimics the physical topology. At machine f, the service rate information kD is passed downstream to the successor while the Kanban arrival rate information cD is passed upstream to all the predecessors. For example, the model at machine 2 receives k information from the machine-4 model, k information from the machine-5
Fig. 7. In#uence diagram for the assembly network in Fig. 1. Note: for f3¸"+7, 8, 9, 10, 11, 12,, kD are already known exactly while cD are not needed.
1038
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
model, and c information from the machine-1 model. The model at machine 2 then generates c information, which it sends to the machine-4 model, c information, which it sends to the machine-5 model, and k information, which it sends to the machine-1 model. Each machine model interacts with all its immediate neighbors, and with no one else, and these information #ows are always bidirectional. The resulting equations are block triangular. For example, the Kanban arrival information c passed to machine 2 contains all the information needed to compute the behavior at the subtree consisting of machines 2 and 4}10. In particular k and k will be implicit functions of c, so that the entire subtree can be collapsed into a single ,ctitious machine 2I. This provides the justi"cation for an aggregation approach whereby the subtrees are merged sequentially, similar to the merging in [7]. The critical feature in any such merging scheme, of course, is making accurate estimates of the queue-dependent service rate kI of the "ctitious machine 2I. Hazra's results show that this approach is highly promising.
7. Summary The performance evaluation of a closed Kanban-controlled assembly system is accomplished by a new heuristic version of the exact aggregation}disaggregation theory for "nite continuous-time Markov chains. It provides accurate estimates of the throughput. Because the approximation is theory based, it provides a framework for further model development, with some possible extensions described in the paper. The approximation has the novel feature of doing simultaneous multiple partitions of the state space, in such a way that the associated aggregate transition rates are mutually consistent. The methodology is a novel approach towards extending aggregation ideas to fork-and-join queueing networks. It provides fast and accurate throughput estimates for complex production systems, and insights into the behavior of these systems. We use the results of this methodology to show how the throughput of an assembly system varies with changes in the parts #ows, in the total number of Kanbans, or in the processing rate at the individual stations. The speed at which we obtain these results make it a practical tool for the analysis of key design trade-o!s for complex manufacturing systems. Anticipated use of the model includes balancing of routes, speeds and topologies.
Acknowledgements Parts of this paper were completed while P. Schweitzer was a Lady Davis Fellow at Hebrew University in Jerusalem, Israel. He thanks the University and Lady Davis Trust for many courtesies extended.
References [1] Hazra J, Seidmann A. Performance evaluation of closed tree-structured assembly systems. IIE Transactions 1996; 28:591}9.
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1039
[2] Takahashi Y. A lumping method for numerical calculations of stationary distributions of Markov chains. Research Report No. B-18, Dept. of Information Sciences, Tokyo Institute of Technology, Tokyo, Japan, 1975. [3] Schweitzer PJ. Aggregation methods for large Markov chains. In: Iazeolla G, Courtois PJ, Hordijk A. editors. Mathematical computer performance and reliability. Amsterdam: North-Holland, 1984:275}86. [4] Schweitzer PJ. A Survey of aggregation}disaggregation in large Markov chains. In: Stewart WJ, editor. Numerical solution of Markov chains, New York: Marcel Dekker, 1991:63}88. [5] Dallery Y, Liu Z, Towsley D, Equivalence, reversibility, Symmetry and concavity properties in fork-join queueing networks with blocking. Journal of Association for Computing Machinery 1994;41:903}42. [6] Viswanadham N, Narahari Y. Performance modeling of automated manufacturing systems. Englewood Cli!s, NJ; Prentice-Hall, 1992. Hoop MJ, Spearman ML. Factory physics, foundation of manufacturing management. Homewood, IL: Irwin, 1996. [7] Hazra J. Performance evaluation of Kanban-controlled assembly networks. PhD Thesis, William E. Simon Graduate School of Business Administration, Rochester, NY: University of Rochester, 1995. [8] Gershwin SB, Manufacturing systems engineering. Englewood Cli!s, NJ: Prentice-Hall, 1994, Buzacott JA, Shanthikumar JG. Stochastic models of manufacturing systems. S. Englewood Cli!s, NJ: Nahmias, Prentice-Hall, 1993 Production and operations analysis. Homewood, IL: Irwin, 1997. [9] Chandrasekhar Rao P, Suri R. Performance analysis of an assembly station with input from multiple fabrication lines. Working paper Dec. 1996, Department of Industrial Engineering, University of Wisconsin-Madison, Madison, WI. [10] Chandrasekhar Rao P, Suri R. Approximate queueing network models for closed fabrication!assembly systems } Part I: Single level systems. Production and Operations Management 1998;3:244}75. [11] Baynat B, Bruno, Dallery Y. Approximate techniques for general closed queueing networks with subnetworks. EJOR 1992; 69(2): 250}64. [12] Schweitzer PJ. Aggregative modelling of queueing networks. In: Bonatti M, editor. Teletra$c science for new cost-e!ective systems, networks and services. Amsterdam: Elsevier, North-Holland, 1989:1522}8. [13] Takahashi Y. Aggregate approximation for acyclic queueing networks with communication blocking. In: Perros HG, Altiok T, editors. Queueing networks with blocking. Amsterdam: Elsevier, 1989:33}46. [14] Dimascolo M, Frein Y, Dallery Y. An analytical method for performance evaluation of Kanban controlled production systems. Operations Research 1996;44(1):50}65.
Jishnu Hazra received his B.Tech (Mech. Engg.) from IIT, Delhi and M.S. and Ph.D. from the University of Rochester in Operations Management. His research interests are in Supply Chain Management and JIT Manufacturing. He has conducted Executive Development Programs on Supply Chain Management for TATA}IBM executives and other companies. He is also involved in consulting and research projects in the areas of R&D Project Planning and Logistics. Paul J. Schweitzer is Professor of Business Administration at the William E. Simon Graduate School of Business Administration of the University of Rochester. His interests include military operations research, performance models of manufacturing and computer systems, Markovian decision processes (MDPs), management science and information technology. He has done seminal work on value-iteration for MDPs, perturbation of Markov chains, mean value analysis of queueing networks including #exible manufacturing systems, and aggregation techniques. He is the author of over 100 research papers and reports. Abraham Seidmann is the Xerox Professor and Areas Coordinator of Computers and Information Systems, Management Science and Operations Management at the William E. Simon Graduate School of Business Administration at the University of Rochester. Professor Seidmann is the author of over 75 research articles which appear in many of the leading scienti"c journals and is a Department Editor on Interdisciplinary Management Research and Applications in Management Science. His current research and consulting activities include Information Systems, Business Process Reengineering, Project Management and Optimal Resource Allocation, Strategic Manufacturing Systems, Information Economics, Healthcare Management, Stochastic Processes and Performance Modeling for Capacity Planning of Assembly Systems.
Analyzing closed Kanban-controlled assembly systems by iterative aggregation}disaggregation Jishnu Hazra , Paul J. Schweitzer *, Abraham Seidmann Indian Institute of Management, Bangalore, 560076, India W.E. Simon Graduate School of Business Administration, University of Rochester, Rochester, NY 14627, USA
Abstract In recent years, many companies have begun using WIP-limiting strategies such as CONWIP, Kanban or drum-bu!er-rope to control parts #ows in complex production systems. This paper analyzes assembly systems with a tree structure, random processing times and a constant WIP control system. We present heuristic version of the exact aggregation}disaggregation theory for "nite Markov chains to evaluate the performance of these closed Kanban-controlled assembly systems. Because the approximation is theory based, it provides a framework for further model development, with some possible extensions described in the paper. The approximation has the novel feature of doing simultaneous multiple partitions of the state space, in such a way that the associated aggregate transition rates are mutually consistent. The methodology is a novel approach towards extending aggregation ideas to fork-and-join queueing networks, and it provides several useful operational and analytical insights. Numerical comparisons with simulation show that the proposed approximation computes accurate estimates of the plant throughput. It provides a fast way to assess the performance and economic impact of changes in the total WIP level (or the number of Kanbans), or in the part routes, on the throughput rate of the assembly system. 1999 Elsevier Science Ltd. All rights reserved. Scope and purpose We develop a new methodology for rapid performance analysis of assembly systems with a tree structure, random processing times, and a constant work in process Kanban control system. Extensive evaluations show that this methodology computes accurate estimates of the plant throughput. It provides a fast way for managers to assess the performance and economic impact of changes in the total WIP level (or the number of Kanbans) and in the part routes, on the overall throughput rate of the assembly system. Keywords: Assembly systems; Kanban-controlled; Aggregation; Disaggregation
* Corresponding author. E-mail address: [email protected] (P.J. Schweitzer) 0305-0548/99/$ - see front matter 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 3 0 5 - 0 5 4 8 ( 9 9 ) 0 0 0 2 4 - 6
1016
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1. Introduction Many manufacturing systems require more than one distinct input part. These machines cannot start processing until all input parts from the preceding stages are available at their input bu!er. This results in synchronization delays, which make such queueing systems di$cult to analyze. The multi-stage assembly system is represented by a directed tree as discussed in [1]. We use the generic term part for components and subassemblies. The part in a bu!er of a downstream machine is an assembly of parts from the bu!ers of upstream machines. We model the detailed parts #ows in the tree-structured assembly system as a continuous-time Markov chain (CTMC). We begin, therefore, with the aggregation}disaggregation (A/D) approach to solving CTMCs with large state spaces and then show how to use it for the rapid approximation of complex assembly queueing systems. We expect our results to become an important part of the process analysis activities undertaken by assembly system designers and managers. The need for our research arises from the observation that even complex queueing network models cannot capture the part mating delays, and simulation models of such systems require special e!ort due to the particular part #ows involved. Our results show that the proposed approximation computes accurate estimates of the assembly plant throughput. 1.1. Exact aggregation}disaggregation We consider an ergodic continuous-time Markov chain with "nite state space S and transition rate matrix +r " i, j3S, iOj,. The equilibrium probabilities [n ] (where n "Pr[in state i]), are GH G GZ1 G the solution to the "S" Kolmogoro! equations: n r " nr ; G GH H HG HZ1 G HZ1 G
i3S (discard one),
(1.1a)
n "1. (1.1b) G GZ1 Often "S" is too large to permit direct solution of Eq. (1.1) either by an exact (pivoting) algorithm or even by successive approximations. An alternative numerical approach is Takahashi's aggregation}disaggregation (A/D) procedure [2}4], in which the state space S is partitioned into J groups of `similara states: S"S(1)#S(2)#2#S(a)#2#S(J).
(1.2)
Attention is focused upon the aggregate equilibrium probabilities, n "Pr[in group S(a)]" n , 1)a)J ? G GZ1? and upon the aggregate transition rates,
n G r , 1)a, b)J, aOb R " GH ?@ n IZ1? I HZ1@ GZ1?
(1.3)
(1.4)
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1017
with the interpretation R dt#o(dt)"Pr[in S(b) at time t#dt" in S(a) at time t], aOb. ?@ [R ] inherits the non-negativity and ergodicity properties from [r ]. ?@ GH The A/D approach replaces Eq. (1.1) with the equivalent quadruple set of equations consisting of Eq. (1.4); the aggregation equations: ( ( na R " n R , 1)a)J (discard one) ?@ @ @? @ @ @$?
(1.5a)
@$?
( n "1; ? ? the disaggregation equations for group S(a), 1)a)J: ( n r " n r # n D G, G GH H HG @ @ @ HZ1 G HZ1? G @$? and the disaggregation rate matrix
n H D G" r , @ n HG IZ1@ I HZ1@ with the interpretation
i3S(a)
(1.5b)
(1.6)
(1.7)
D G dt#o(dt)"Pr[in i at time t#dt" in S(b) at time t]. @ Formally, Eq. (1.5) is the set of Kolmogoro! equations for an aggregated J-state CTMC, even though the aggregated process is not Markovian. The quadruple of Eqs. (1.4)}(1.7) is just su$cient to determine the quadruple of unknowns +n , n , R , D G,. The usual iterative A/D algorithm G ? ?@ @ proceeds by alternation between solving Eq. (1.5) for +n , and solving Eq. (1.6) for +n ,. Much ? G smaller sets of linear equations are solved at each step. Upon convergence, Eq. (1.3) provides a numerical check. The A/D framework, however, su!ers from two drawbacks. First, since the size of the problem now is even larger than "S", it is still impractical for problems with enormous state spaces. Second, if the system performance characteristics, such as throughputs, depend only on +n ,, it seems ? cumbersome to have to deal with both +n , and +n , in order to compute just the latter. G ? Ideally one wants a purely aggregated version of the problem, one which involves just the +n , ? and the +R ,. This would overcome both drawbacks. But a glance at the quadruple shows that this ?@ is impossible, because the +n , depend upon the +R ,, which in turn depend on the +n ,. The exact ? ?@ G +n , cannot be computed, therefore, without already knowing +n ,. ? G 1.2. Heuristic aggregation The pursuit of a purely aggregated version, involving just +n , and +R ,, may however, be ? ?@ undertaken as an approximation. One estimates +R , and then solves Eq. (1.5) to obtain an estimate ?@ of +n ,. As the estimate of +R , improves, so does the estimate of +n ,. ? ?@ ?
1018
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
A whole nested family of procedures can be used to estimate +R ,, with increasing accuracy ?@ accompanied by increasing computational e!ort. At one extreme is solving the full quadruple set of equations with perfect accuracy. At the other extreme is a one-pass procedure, in which a single estimate of the conditional probabilities [n / n ] is inserted into Eq. (1.4) to obtain a G IZ1? I single estimate of the +R ,. In between are various heuristic estimates of +R , that may ?@ ?@ involve intermediate models, typically iterative and more coarsely aggregated. These heuristics exploit structural features of the problem at hand, and no general recipe can be given for their construction. 1.3. Multiple partitions This paper describes one such intermediate procedure as applied to closed Kanban-controlled assembly systems. It has the novel feature of doing multiple partitions of the state space, and it generates estimates of +R , by the requirement that the various partitions generate mutually ?@ consistent estimates of the aggregate transition rates. This leads to a xxed-point problem for the +R ,, which itself is solved by iteration. We show that the procedure works, and discuss how it is ?@ capable of even higher accuracy by additional re"nement of the procedure for estimating the +R ,. ?@ The plan of the paper is as follows. Section 2 describes the Kanban system and our notation. Section 3 sets up a single partition of the state space S. Section 4 sets up the multiple partitions and the resulting "xed-point problem. Section 5 describes the computational procedure and results, and Section 6 o!ers further development of the model. Section 7 contains our summary and conclusions.
2. Closed Kanban-controlled assembly system In this paper we study an assembly system operating under CONWIP (Constant Work In Process) control. The goal is to estimate the throughput and the expected local bu!er occupancy for a given topology and a given number of Kanbans. The Kanban system dynamics can be understood from Fig. 1. The assembly system is a directed tree with M machines numbered 1, 2,2, M (here, M"12). Machine f has exponentially distributed assembly times with mean assembly time 1/k , 1)f)M. There is one root node r (here, r"1), a set ¸ of two or more leaf D nodes (here, ¸"+7, 8, 9, 10, 11, 12,), and a set INT of intermediary nodes (here, INT"+2, 3, 4, 5, 6,). New components are loaded at the leaf nodes, processed or assembled at each machine, and forwarded to the machine's output bu!er. For instance, machine 8 takes parts from the input bu!er (whose queue length is n ) and, after completion, forwards them to its output bu!er (whose queue length is n ). That is, n will decrease by one, and n will increase by one. Leaf nodes start the production process of the basic platform (e.g., a computer's motherboard or a car's chassis.) These nodes do not do assembly. Note that machine 5 takes one part from each of its two input bu!ers (whose queue lengths are n and n ) and after assembly forwards them to its output bu!er (whose queue length is n ). Both n and n will decrease by one, and n will increase by one. The root node, machine 1, combines subassemblies coming from nodes 2 and 3. When the root node completes assembly, a fully-assembled part exits the system, and a Kanban (authorization to
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1019
Fig. 1. Part #ows and bu!er structure in a 12-machine assembly system.
produce) as well as raw material is instantaneously and simultaneously delivered to all the leaf nodes. That is, n and n will both decrease by one, while n , n , n , n ,n , and will all increase by one. n In this notation, n denotes the number of parts in the bu!er located downstream from i and GH upstream from j. n includes both parts on queue and parts in the midst of service at j. There is GH su$cient bu!er space that blocking never occurs. Starvation, however, can occur; for example, machine 5 cannot begin assembly unless it has at least one part in each of its input bu!ers, namely, when n *1 and n *1. Aside from the feedback #ow of Kanbans out of the root node, the topology is a directed tree. In general, the tree is unbalanced. Every machine has at least one input bu!er and, aside from the root node, exactly one output bu!er. Technically, this is a closed queueing network with one fork and several joins [5]. The Kanban control ensures that balanced and coordinated #ows of raw materials occur at all the leaf nodes. In open systems, this is often accomplished by MRP systems. A special case of the tree system shown here is the linear assembly line, with one leaf node, which has been analyzed as a CONWIP (constant work-in-process) system. For extensive discussion of models of closed Kanban systems, including CONWIP, see [1,6,7]. Other recent discussions are given in [5,8}11,14]. There is an extensive literature on open assembly systems, but not on closed ones. A convenient modeling viewpoint is that parts never disappear; they are merely mated together. For example a Kanban will simultaneously bring 6 components to the 6 leaf nodes. Four of these six components will subsequently enter the root node from machine 2, while two others will subsequently enter the root node from machine 3. The completed part consists of the 6 assembled components. After exiting the root node, the completed part is `disassembleda into its 6 components, and these 6 components are delivered back to the 6 leaf nodes. The total number of components never changes.
1020
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
The assembly system is started with all bu!ers empty except for the leaf nodes, each of which have their input bu!ers "lled with N components. There will then always be enough components to assemble exactly N "nished goods. We let m denote the number of parts going upstream on a path from machine f to a leaf node. D A convenient consequence of the above assumptions is path independence; m is independent of the D choice of path between f and a leaf node. In Fig. 1, this means: m "n #n "n #n (equality of two paths) m "n #n #n #n "n #n #n #n "n #n #n "etc. (equality of 6 paths)
(2.1)
"N (Note that m "N). P A convenient way of writing the property of path independence is as the recursion m "m #n fOr, where S( f ) denotes the successor to machine f. 1D D D1D The boundary values are m "n if f is a leaf. For Fig. 1, by taking f"8 or 9, D PD m "m #n "m #n . Thus the two paths between machine 5 and the two leaves 8 and 9 have the same number of parts. Path independence is lost if any machine requires more than one part from any of its input bu!ers, or if initial conditions have unequal numbers of parts on the di!erent paths. For simplicity our analysis employs these convenient conditions, but it could be undertaken without them. The network topology is characterized by M"number of nodes, r"root node, INT"set of intermediate nodes, neither root nor leaf, ¸"set of leaf nodes or nodes without predecessors (leaf nodes have only one input bu!er), S( f )"unique successor (downstream) node to node f, fOr [e.g., S(8)"5], P( f )"non-empty set of predecessor (upstream) nodes to node f, f,¸ [e.g., P(1)"+2, 3,]. This topological information, plus N and the +k ,, constitute the full set of parameters for the G problem. The full state space S consists of the set of bu!er occupancy numbers n"+n " 1)i)M, j3S(i),, restricted by the equality of parts on each path and by having N parts GH on each leaf-to-root path. The dynamics can be described as a continuous-time Markov chain (CTMC) with states n"+n , and state space S. The transition rate for machine i (i3INT) to complete an assembly is GH k if n *1 for all j3P(i), G HG 0 otherwise.
If such a transition occurs, the output bu!er level n S(i) increases by one while all input bu!er levels G +n " j3P(i), decrease by one. HG
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1021
The equilibrium probabilities of this CTMC are denoted by n(n), n3S. Since realistic problems can easily have 10 states, it is hopeless to try to calculate them directly from Eq. (1.1). Furthermore, they are not of central interest. The main performance measures of interest are h "throughput of node f, 1)f)M, D calculated from h "k Pr [all input bu!ers for node f are non-empty;] D D U "utilization of node f, 1)f)M, D calculated from
(2.2)
; "h /k ; (2.3) D D D E[n ]"mean occupancy level at the output buwer of node i; 1)i)M, iOr, j"S(i) GH E[n ]"mean occupancy level at the input buwer to each leaf node j, j3L. PH These are calculated from E[n ]" n n(n). GH GH n
(2.4)
The problem can be made dimensionless by measuring the service rates +k , and the throughputs D +h ,, as multiples of k or of min k . D P E E The following properties, Eqs. (2.5) and (2.6), will hold for the exact model and should hold, at least as approximate equalities, for any approximation: h "h "2"h . (2.5) + All machines must have the same throughput, and this cannot exceed min k . Due to WDW+ D starvation, the throughput can be appreciably less than min k . D D The path-independence condition of mean bu!er levels upstream of any machine [e.g., Eq. (2.1)]. (2.6) A more systematic way of writing the path-independence condition is as the recursion Em "Em #En S( f ), all fOr 1D D D plus the boundary values Em "N and Em "En (if f is a leaf). P D PD
(2.7)
3. The partition associated with machine f 3.1. The aggregate states Pick any non-leaf machine f with two or more input bu!ers. For simplicity the description is given as in Fig. 2 for the case in which f has 2 input bu!ers and two predecessors, labeled g and h.
1022
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
Fig. 2. Bu!er content feeding an assembly mode.
[The analysis is simpler if f has one input bu!er, and more cumbersome if f has 3 or more input bu!ers. The framework is the same in all cases.] We generate a partition based on the triple (m , n , n ), of non-negative integers, where D ED FD 0)max[n , n ])m )N. ED FD D That is, a in Eq. (1.2) is this triple a"(m , n , n ), D ED FD and S(a) is the set of states n3S with the given values of m , n , n . All other state information is D ED FD discarded. The aggregate probability [see Eq. (1.3)] that this set of states occurs is denoted by nD(m , n , n ), with the superscript indicating this as a partition based on machine f. S(a) is D ED FD distinguished from S( f ) because a is a triple while f is an integer between 1 and M. It is clear that any accurate analysis of machine f must contain at least information about both of its input bu!er levels n and n . Carrying some additional upstream information via m is useful ED FD D because it supplies m "m !n and hence some partial information about how likely g is to E D ED complete a service, thereby increasing n by one. ED 3.2. Aggregate transition rates According to Eq. (1.5), for the purpose of computing aggregate transition rates, the changes in the triple a"(m , n , n ) behave as a CTMC. The four possible transitions for the triple are as D ED FD follows for f3INT: (a) Machine f has a service completion Next triple b:m , n , n all decrease by one. The aggregate transition rate `R a is D ED FD ?@ k I (n *1 and n *1) and is known. Here I(E) is the indicator function, 1 if E holds and D ED FD 0 otherwise. (b) Machine g has a service completion Next triple b:n increases by one, m and n are unchanged. ED D FD Case 1: g is a leaf. The aggregate transition rate `R a is k I (m !n *1) and is known. ?@ E D ED Case 2: g is not a leaf. The aggregate rate `R a is unknown and will have to be estimated. ?@
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1023
Denote the exact value by kE(m , n , n ). The expression kE(m , n , n ) is to be interpreted D ED FD D ED FD as the e!ective service rate of g as seen by f when in the aggregate state a"(m , n , n ). The D ED FD rate is given by Eq. (1.4) as kE(m , n , n )"k Pr [all of g s input bu!ers are non-empty " m , n , n ] D ED FD E D ED FD In both cases, the rate is zero unless m "m !n *1. E D ED (c) Machine h has a service completion Analogous to (b): if h is not a leaf, there is an unknown e!ective service rate kF(m , n , n ) for D ED FD transition to a new state where n increases by one. FD (d) Kanban arrives simultaneously at all the leaf nodes The next triple b has m increased by one, with n and n unchanged. The aggregate rate D ED FD `R a is unknown and must be estimated. We denote the exact aggregate transition rate `R a ?@ ?@ as cD(m , n , n ). It is to be interpreted as the arrival rate of Kanbans as seen by f while in D ED FD aggregate state a"(m , n , n ). D ED FD The rate is given by Eq. (1.4) as cD(m , n , n )"k Pr[all input bu!ers to root node are non-empty " m , n , n ]. D ED FD P D ED FD Of course this rate vanishes if fOr and m "N, because the root node is starved for parts D coming from machine f. If f3INT, the Kanban arrival rate cD(m , n , n ) is unknown and must be estimated. If f"r, D ED FD the Kanban arrival rate is known exactly: cP(m "N, n , n )"k I(n *1 and n *1). P EP FP P EP FP The aggregate CTMC associated with a non-leaf machine f has as its states the triples (m , n , n ) D ED FD with 0)max[n , n ])m )N, and it involves some unknown aggregate transition rates ED FD D kE, kF and cD. If these rates were known exactly, solution of this aggregate problem (1.5) would give exact values for the aggregate equilibrium probabilities: nD(m , n , n )"Pr[in a state with given values of m , n , n ]. D ED FD D ED FD 3.3. Approximate aggregate transition rates The aggregate transition rates are unknown and must be estimated. We denote estimates with a*, namely kE*, cD*, nD*, h*, Pr[ ]*, En* , and En* . We use approximations of the form (3.1) and D ED FD (3.5) below. Speci"cs about kE* and cD* are given in Section 4.1. The key approximation for kE is kE(m , n , n ) kE*(m !n ); g3INT D ED D ED FD
(3.1)
i.e., the e!ective service rate at g depends only on m . If g is a leaf node, Eq. (3.1) is replaced by the E exact result: kE(m , n , n )"k I(m !n *1), g3¸. D ED FD E D ED
1024
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
We know kE(m , n , n )"0 if m "0, and we therefore require that D ED FD E kE*(m )"0 if m "0. E E
(3.2)
In addition, since kE(m , n , n )"k Pr[g busy " m , n , n ], D ED FD E D ED FD
(3.3)
we know that kE(m , n , n ))k , and we therefore require that D ED FD E kE*(m ))k , m "1, 2,2, N. E E E
(3.4)
The key approximation for cD is cD(m , D cD(m , D
n , ED n , ED
n ) cD*(m ) fOr, D FD n )"k I(n *1 and n *1) f"r, FD P EP FP
(3.5)
i.e., the e!ective Kanban arrival rate seen by fOr depends only on m . D Since cD(n , n , n )"0 when fOr and m "N (the root machine is starved for parts from D ED FD D machine f ), we therefore require that cD*(N)"0, fOr.
(3.6)
cD(m , n , n )"k Pr[r busy " m , n , n ] D ED FD P D ED FD
(3.7)
Since
(The Kanban arrival rate " the completion rate at the root), we know that cD(m , n , n ))k , D ED FD P and we therefore require that cD*(m ))k . D P
(3.8)
3.4. Special cases for f If machine f has only one input bu!er, then the partition is on the pair (m , n ). The analysis is D ED simpler because transitions of type (c) are absent. If f is the root machine r, the analysis is somewhat simpler because m "N, so only n and P EP n enter the CTMC. In this case, transitions (a) and (d) coincide and have rate FP k I(n *1 and n *1) to a new triple where both n and n have decreased by one. P EP FP EP FP 3.5. Model outputs Given these estimated transition rates +kE*, kF*, cD*,, we can solve the aggregate CTMC (1.5) and obtain estimates nD*(m , n , n ) of the exact aggregate equilibrium probabilities D ED FD nD(m , n , n )["`naa]. D ED FD
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1025
3.6. Byproducts of partition at f As by products, we can calculate estimates of throughputs and mean queue lengths from Eqs. (3.9) and (3.10): h*"estimated throughput at f D "k Pr[n *1 and n *1]* [cf. Eq. (2.2)] D ED FD , KD KD "k nD*(m , n , n ) f3¸. (3.9a) D D ED FD KD LED LFD If f is a leaf node, there is no aggregate model such as in Fig. 2 to use for computing h , as in Eq. D (3.9a). However, f is part of the aggregate model
¹"set of predecessors of S( f ) other than f, i.e. siblings of f. associated with S( f ) and we may use h*"h* , f3¸. D 1D En* "estimated queue length at output bu!er of g ED , KD KD " n nD*(m , n , n ). ED D ED FD KD LED LFD
(3.9b)
(3.10)
4. Mutually consistent rates 4.1. Estimating the aggregate rates At this point we have constructed a partition and an aggregate model for one non-leaf machine f, with predecessors g and h, with model inputs kE*, kF* and cD*, and with model outputs nD*(n , n , n ), h*, En* , and En* . Now consider performing a similar partition on every D ED FD D ED FD non-leaf machine. We must provide the model}inputs in a self-consistent way. This will lead to a joint "xed-point problem [Eqs. (4.1)}(4.3) below] for the full set of unknown rates, +kD*(m ), 1)m )N,#+cD*(m ), 0)m )N!1, for all f3INT. D D D D The number of unknowns is 2N"INT""O(NM), and the storage requirements are modest for realistic problems. 4.2. Estimating the ewective service rates The e!ective service rates are estimated as kD*(m )"service completion rate at machine f, given that there are m jobs upstream of f D D "k Pr[n *1 and n *1"m ]*, D ED FD D
1026
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
so kD*(m )"0, D
m "0 D
k KD KD nD*(m , D kD*(m )" D LED LFD D KEDD KFDD nD*(m , L L D
n , ED n , ED
n ) FD , n ) FD
1)m )N D
f3INT, P( f )"+g, h,.
(4.1)
This automatically satis"es Eqs. (3.2) and (3.4).
4.3. Estimating the Kanban arrival rate The Kanban arrival rate, cE*(m ), is estimated by the consistency requirement that g and its E successor f"S(g) `seea the same expected Kanban arrival rate, under the same conditions, namely the same m : E cE*(m )"E[cD*(m ) " m ], 0)m )N, f"S(g) E D E E
0
"
,D EcD (m )Pr[m " m ]* D D E K K *
m "N E fOr, g3INT, 0)m )N!1. E
The last expression is cE*(m )"0, m "N, E E cD*(m )nD*(m , n , n ) D D ED FD , 0)m )N!1, where f"S(g)Or. cE*(m )" E E nD*(m , n , n ) D ED FD
(4.2)
The triple sums are taken over all triples (m , n , n ) which satisfy m )m )N, 0)n )m , D ED FD E D ED D 0)n )m , m !n "m . FD D D ED E When f"r and g3P(r), Eq. (4.2) is replaced by cE*(m )"k Pr[machine r is busy " m ]* E P E "k Pr[n *1 and n *1 " m ]* P EP FP E
0
" k ,\KE nP (m "N, n "N!m , n "k) P I P EP E FP ,\KE nP*(m "N, n "N!m , n "k) I P EP E FP *
m "N, E
(4.3)
0)m )N!1. E
Note that Eq. (4.3) used the exact expression c P(m , n , n )"k I(n *1 and n *1) to evaluate P EP FP P EP FP cE*. Eqs. (4.2) and (4.3) automatically satisfy Eq. (3.6). Eq. (4.3) satis"es (3.8). Induction on the nodes, proceeding upstream from the root, then shows that Eq. (4.2) satis"es (3.8).
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1027
4.4. Overall xxed-point problem Eqs. (4.1)}(4.3) for all f3INT#+r, provide a set of simultaneous relationships among the unknowns +kD*(m ), 1)m )N, f3INT,#+cD*(m ), 0)m )N!1, f3INT,. D D D D They constitute a xxed-point problem for these unknowns. As is the case for other nodal decomposition schemes, this "xed-point set can be solved by successive approximations: kD*(m )"k I(m *1), f3INT#+r,, 0)m )N, D D D D cD*(m )"k I(m (N), f3INT#+r,, 0)m )N. D P D D Loop step: pick next machine f3INT#+r, and generate partition associated with f retrieve the aggregate rates cD* and +kE*, g3P( f ), solve aggregate CTMC for nD* revise kD* [if f3 INT] by Eq. (4.1) revise cE* [if g3P(f) such that g,¸] by Eq. (4.2) or Eq. (4.3) revise h* by Eq. (3.9) D revise En * for all g3P( f ) by Eq. (3.10) ED Termination test: if revised values have not converged, return to loop step. Initialization:
This algorithm typically converges in a few score iterations. The most time-consuming part is the repeated solution of the aggregate CTMCs. We did not explore acceleration by exploiting the skip-free behavior of m , i.e., a generalized birth and death D process. The aggregate CTMCs have O(NR>) states where t"max "P( f )". Even with small t, this D puts a signi"cant restriction on N. When one moves onward from the current choice of f to the next, the nD* can be discarded. Only the +kD*, and +cD*, need be carried forward, so the algorithm's storage requirements are only O(MN). Storing the nD* has the merit of providing a good starting point for the iterative solution of this aggregate CTMC when one returns to it, however, thus reducing computer time. The scheme merits the name iterative aggregation}disaggregation because it alternates between computing the aggregate parameters kD* and cD*, and solving the detailed Kolmogoro! equations for one partition, +nD*(m , n , n ),. Strictly speaking, both of these steps are at the D ED FD aggregate level, but the former is more highly aggregated than the latter. Upon termination, we immediately have the performance measures +h*, and +En* ,. Other D ED performance measures can be computed from nD*.
5. Computational results We start by presenting our results for a 7-machine assembly system with binary tree topology (Fig. 3). It was evaluated both by simulation and by the aggregation}disaggregation algorithm for populations N"6, 9, 12. Table 1 displays the initial numerical results. The principal observations are that the algorithm successfully computes throughput to within 1}2%. It has much lower accuracy (up to 30% error) for predicting expected queue lengths. From a practical system-design perspective, it is important to note that the absolute error in the mean queue length hardly exceeds
1028
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
Fig. 3. The experimental seven-machine assembly system. Every k "1. G
Table 1 Numerical results for test case N
h En "En En "En "En "En En "En "En "En
6
9
12
A/D
Sim
%error
A/D
Sim
%error
A/D
Sim
%error
0.658 2.220 2.600 1.180
0.664 2.360 2.140 1.500
!1.0 !6.0 22.0 !22.0
0.750 3.270 3.900 1.820
0.750 3.490 3.210 2.300
0.0 !3.5 22.0 !20.0
0.804 4.320 5.350 2.430
0.801 4.630 4.270 3.080
0.4 !7.0 23.0 !21.0
Note: h"throughout, En "expected occupancy in output bu!er of machine i. GH %error"100(algorithm!simulation)/simulation. Mean queue lengths are given as average of values at these bu!ers, which should be identical. h"average of the h , which should be identical. D
one part. We attribute this inaccuracy to the inability of the aggregate model to describe correlated behavior among machines. This defect is common to most single-node decomposition methods [12]. This may not be a major practical issue as here we deal with closed queueing systems where bu!er occupancies are also controlled through the size of the population. The assembly system in Fig. 3 has four loops of Kanbans (and parts). We expect, therefore, four relationships of the form En #En #En "N. These held within 1%. The solution violated Eq. (2.5), but the variation among the h* was only about 1%. The average of the h* seems to be D D a good estimate of h. Similarly, the solution violated Eq. (2.7), but the variation among di!erent paths was small as well. Section 6.2 sketches one way to force equalities in (2.5) and (2.7). We "nd that the assembly system throughputs are only 2/3 as large as the ideal value, min +k ,"1, for N"6. That happens because of inherent delays among mating parts as evidD G enced by the increased bu!er occupancies in front of those machines. The results in Table 1 show that the system throughput increases in a concave fashion to 4/5 of the ideal value as we double the number of Kanbans (N"12).
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1029
Table 2 Input data for the seven-machine network in Fig. 1 Data set k
k
k
k
k
k
k
1 2 3 4
1 0.8 0.8 1
1 0.8 1 1
1 0.8 1 0.8
1 0.8 1 0.8
1 0.8 0.8 0.8
1 0.8 0.8 0.8
1 1 1 1
Table 3 Input data for the six-machine network in Fig. 4 Data set k
k
k
k
k
k
1 2 3 4
1 1 1 1
1 0.8 0.8 1
1 1 1 0.8
1 1 1 0.8
1 1 0.8 1
1 1 1 1
Table 4 Input data for the "ve-machine network in Fig. 5 Data set
k
k
k
k
k
1 2 3
1 1 1
1 0.8 1
1 0.8 1
1 1 0.8
1 1 0.8
Table 5 Input data for the "ve-machine network in Fig. 6 Data set
k
k
k
k
k
1 2 3 4
1 1 1 1
1 1 1 1
0.8 0.8 0.8 1
0.8 1 1 1
0.8 0.8 1 1
We conducted extensive computations to better understand the accuracy performance of the algorithm. The results reported here are based on four di!erent network topologies (Figs. 1 and 4}6). For each topology, three or four di!erent data sets are used with three di!erent values of N, the number of Kanbans (N"6, 9 and 12). Kanban #ows from the root node (always numbered 1 in the diagram) to the leaf nodes. The data sets are given by Tables 2}5. Our results in Tables 6}19
1030
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
Table 6 Results for the seven-machine network and Data set 2 in Table 2 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.540 1.977 1.977 2.751 2.751 2.751 2.751 1.275 1.275 1.275 1.275
0.558 2.048 1.968 2.347 2.258 2.338 2.386 1.606 1.694 1.695 1.647
!3.23 !3.47 0.46 17.21 21.83 17.66 15.30 !20.61 !24.73 !24.78 !22.59
0.613 2.803 2.803 4.181 4.181 4.181 4.181 2.008 2.008 2.008 2.008
0.626 2.954 2.735 3.570 3.421 3.576 3.721 2.477 2.625 2.689 2.544
!2.08 !5.11 2.49 17.11 22.22 16.92 12.36 !18.93 !23.50 !25.33 !21.07
0.651 3.602 3.602 5.661 5.661 5.661 5.661 2.727 2.727 2.727 2.727
0.666 3.807 3.501 4.850 4.595 4.860 5.050 3.344 3.598 3.639 3.449
!2.25 !5.38 2.88 16.72 23.20 16.48 12.10 !18.45 !24.21 !25.06 !20.93
Table 7 Results for the seven-machine network and Data set 3 in Table 2 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.587 2.235 1.934 2.804 2.804 2.544 2.544 0.962 0.962 1.523 1.523
0.600 2.373 1.941 2.445 2.389 2.023 2.100 1.182 1.238 2.037 1.959
!2.17 !5.82 !0.36 14.68 17.37 25.75 21.14 !18.61 !22.29 !25.23 !22.26
0.662 3.239 2.760 4.363 4.363 3.753 3.753 1.402 1.402 2.489 2.489
0.670 3.544 2.713 3.808 3.727 2.946 3.103 1.648 1.729 3.341 3.185
!1.19 !8.61 1.73 14.57 17.06 27.39 20.95 !14.93 !18.91 !25.50 !21.85
0.707 4.108 3.628 6.091 6.091 4.849 4.849 1.800 1.800 3.526 3.526
0.708 4.735 3.409 5.218 5.149 3.822 4.079 2.047 2.117 4.768 4.512
!0.14 !13.24 6.42 16.73 18.29 26.87 18.88 !12.07 !14.97 !26.05 !21.85
show that the algorithm computes the system throughput with an error of 5% or less in all cases. That error was not a!ected by the value of N. Most of the estimates for the bu!er occupancy were accurate within 30%, or better and within roughly 1 job in absolute value. Estimation errors for bu!er occupancies were greater for upstream machines than for downstream machines. Again, the magnitude of that estimation error was not necessarily correlated with N. In all cases, increasing N results in a concave increase in the system throughput h. The e!ect of changes in N on h depends of the topology of the assembly system, as can be seen by comparing the two "ve-machine topologies studied. The results of our comparative study are presented in Tables 6}19.
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1031
Table 8 Results for the seven-machine network and Data set 4 in Table 2 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.575 2.063 2.063 2.479 2.479 2.479 2.479 1.460 1.460 1.460 1.460
0.592 2.133 2.074 2.005 1.914 1.960 2.031 1.862 1.953 1.966 1.895
!2.87 !3.28 !0.53 23.64 29.52 26.48 22.06 !21.59 !25.24 !25.74 !22.96
0.650 2.956 2.956 3.656 3.656 3.656 3.656 2.391 2.391 2.391 2.391
0.663 3.072 2.919 2.953 2.741 2.861 3.021 2.975 3.187 3.220 3.061
!1.96 !3.78 1.27 23.81 33.38 27.79 21.02 !19.63 !24.98 !25.75 !21.89
0.694 3.808 3.808 4.799 4.799 4.799 4.799 3.397 3.397 3.397 3.397
0.703 3.925 3.707 3.829 3.552 3.698 3.995 4.246 4.523 4.596 4.298
!1.28 !2.98 2.72 25.33 35.11 29.77 20.13 !20.00 !24.89 !26.09 !20.96
Table 9 Results for the six-machine network of Fig. 4 and Data set 1 in Table 3 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.691 1.910 2.539 2.784 2.784 1.803 1.304 1.304 1.642
0.683 2.216 2.751 2.235 2.175 1.639 1.549 1.609 1.611
1.17 !13.81 !7.71 24.56 28.00 10.01 !15.82 !18.96 1.92
0.775 2.842 3.664 4.161 4.161 2.719 1.996 1.996 2.604
0.766 3.309 4.040 3.357 3.263 2.489 2.335 2.428 2.470
1.17 !14.11 !9.31 23.95 27.52 9.24 !14.52 !17.79 5.43
0.825 3.822 4.699 5.552 5.552 3.636 2.624 2.624 3.656
0.813 4.414 5.371 4.453 4.353 3.319 3.134 3.234 3.311
1.48 !13.41 !12.51 24.68 27.54 9.55 !16.27 !18.86 10.42
6. Further re5nements We sketch several ways that the accuracy of the model can be improved with increased computational e!ort. This demonstrates the #exibility of the A/D framework for developing whole families of approximations. 6.1. More-detailed partitions The current partitioning (m , n , n ) on f and its immediate predecessors, g and h, does not D ED FD provide information either upstream of g and h or downstream of f. A "ner partition would employ
1032
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
Table 10 Results for the six-machine network of Fig. 4 and Data set 2 in Table 3 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.654 2.335 2.006 2.535 2.535 2.498 1.130 1.130 1.483
0.647 2.579 2.154 2.022 1.969 2.401 1.399 1.452 1.445
1.08 !9.46 !6.87 25.37 28.75 4.04 !19.23 !22.18 2.63
0.728 3.695 2.650 3.661 3.661 4.121 1.643 1.643 2.220
0.719 4.131 2.822 2.881 2.798 4.102 1.989 2.072 2.076
1.25 !10.55 !6.09 27.07 30.84 0.46 !17.40 !20.70 6.94
0.768 5.195 3.133 4.718 4.718 5.888 2.087 2.087 2.971
0.755 5.942 3.248 3.579 3.490 6.176 2.479 2.568 2.576
1.72 !12.57 !3.54 31.82 35.19 !4.66 !15.81 !18.73 15.33
Table 11 Results for the six-machine network of Fig. 4 and Data set 3 in Table 3 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.623 2.638 1.587 2.352 2.352 2.270 1.009 1.009 2.133
0.615 2.892 1.743 1.844 1.796 2.125 1.265 1.313 2.133
1.30 !8.78 !8.95 27.55 30.96 6.82 !20.24 !23.15 0.00
0.689 4.292 1.956 3.303 3.303 3.506 1.405 1.405 3.534
0.681 4.721 2.202 2.541 2.455 3.389 1.738 1.824 3.409
1.17 !9.09 !11.17 29.99 34.54 3.45 !19.16 !22.97 3.67
0.722 6.139 2.232 4.136 4.136 4.769 1.725 1.725 4.999
0.714 6.823 2.497 3.060 3.003 4.737 2.117 2.174 4.766
1.12 !10.02 !10.61 35.16 37.73 0.68 !18.52 !20.65 4.89
Table 12 Results for the six-machine network of Fig. 4 and Data set 4 in Table 3 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.625 1.372 3.113 2.817 2.817 1.525 1.811 1.811 1.340
0.627 1.631 3.236 2.232 2.134 1.397 2.138 2.236 1.366
!0.32 !15.88 !3.80 26.21 32.01 9.16 !15.29 !19.01 !1.90
0.691 1.799 4.866 4.158 4.158 2.148 3.042 3.042 1.964
0.696 2.158 5.067 3.347 3.105 1.980 3.495 3.737 1.953
!0.72 !16.64 !3.97 24.23 33.91 8.48 !12.96 !18.60 0.56
0.729 2.179 6.713 5.304 5.304 2.748 4.516 4.516 2.522
0.731 2.548 7.217 4.398 4.045 2.418 5.054 5.408 2.366
!0.27 !14.48 !6.98 20.60 31.12 13.65 !10.65 !16.49 6.59
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1033
Table 13 Results for the "ve-machine network of Fig. 5 and Data set 1 in Table 4 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.733 2.189 2.189 1.968 1.968 1.829 1.829
0.701 2.670 2.534 1.679 1.737 1.651 1.730
4.56 !18.01 !13.61 17.21 13.30 10.78 5.72
0.816 3.153 3.153 2.897 2.897 2.938 2.938
0.781 3.943 3.754 2.569 2.620 2.488 2.627
4.48 !20.04 !16.01 12.77 10.57 18.09 11.84
0.867 4.060 4.060 3.753 3.753 4.179 4.179
0.826 5.230 5.017 3.439 3.461 3.331 3.523
4.96 !22.37 !19.08 9.13 8.44 25.46 18.62
Table 14 Results for the "ve-machine network of Fig. 5 and Data set 2 in Table 4 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.650 2.125 2.125 2.424 2.424 1.437 1.437
0.632 2.448 2.309 2.215 2.297 1.337 1.394
2.85 !13.19 !7.97 9.44 5.53 7.48 3.08
0.715 3.052 3.052 3.866 3.866 2.071 2.071
0.697 3.531 3.325 3.653 3.736 1.816 1.939
2.58 !13.57 !8.21 5.83 3.48 14.04 6.81
0.754 3.923 3.923 5.407 5.407 2.662 2.662
0.731 4.542 4.284 5.234 5.348 2.225 2.368
3.15 !13.63 !8.43 3.31 1.10 19.64 12.42
Table 15 Results for the "ve-machine network of Fig. 5 and Data set 3 in Table 4 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.650 2.004 2.004 1.646 1.646 2.337 2.337
0.631 2.453 2.292 1.358 1.410 2.189 2.298
3.01 !18.30 !12.57 21.21 16.74 6.76 1.70
0.710 2.828 2.828 2.254 2.254 3.906 3.906
0.698 3.572 3.224 1.896 1.978 3.533 3.798
1.72 !20.83 !12.28 18.88 13.95 10.56 2.84
0.744 3.577 3.577 2.701 2.701 5.708 5.708
0.732 4.652 4.117 2.329 2.393 5.019 5.491
1.64 !23.11 !13.12 15.97 12.87 13.73 3.95
1034
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
Table 16 Results for the "ve-machine network of Fig. 6 and Data set 1 in Table 5 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.609 1.574 3.613 2.719 2.719 2.387 1.709 1.709
0.624 1.645 3.652 2.244 2.118 2.348 2.112 2.238
!2.40 !4.32 !1.07 21.17 28.38 1.66 !19.08 !23.64
0.672 2.327 5.047 3.935 3.935 3.952 2.740 2.740
0.688 2.332 5.176 3.301 3.011 3.824 3.367 3.657
!2.33 !0.21 !2.49 19.21 30.69 3.35 !18.62 !25.08
0.706 3.184 6.091 4.945 4.945 5.913 3.870 3.870
0.722 3.012 6.482 4.270 3.793 5.518 4.718 5.196
!2.22 5.71 !6.03 15.81 30.37 7.16 !17.97 !25.52
Table 17 Results for the "ve-machine network of Fig. 6 and Data set 2 in Table 5 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.635 1.756 3.377 3.146 2.361 2.623 1.103 1.881
0.645 1.841 3.447 2.777 1.750 2.553 1.383 2.410
!1.55 !4.62 !2.03 13.29 34.91 2.74 !20.25 !21.95
0.700 2.682 4.556 4.748 3.296 4.444 1.575 3.020
0.707 2.739 4.743 4.363 2.301 4.257 1.898 3.961
!0.99 !2.08 !3.94 8.82 43.24 4.39 !17.02 !23.76
0.733 3.738 5.412 6.112 4.322 6.591 2.142 3.936
0.736 3.656 5.799 6.061 2.666 6.201 2.283 5.678
!0.41 2.24 !6.67 0.84 62.12 6.29 !6.18 !30.68
Table 18 Results for the "ve-machine network of Fig. 6 and Data set 3 in Table 5 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.673 2.090 3.001 2.680 2.680 2.999 1.232 1.232
0.679 2.224 3.107 2.238 2.150 2.893 1.538 1.626
!0.88 !6.03 !3.41 19.75 24.65 3.66 !19.90 !24.23
0.742 3.438 3.697 3.810 3.810 5.303 1.746 1.746
0.744 3.659 3.817 3.168 3.015 5.183 2.173 2.326
!0.27 !6.04 !3.14 20.27 26.37 2.32 !19.65 !24.94
0.782 4.810 4.320 4.962 4.962 7.690 2.225 2.225
0.773 5.470 4.107 3.893 3.677 7.893 2.637 2.853
1.16 !12.07 5.19 27.46 34.95 !2.57 !15.62 !22.01
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1035
Table 19 Results for the "ve-machine network of Fig. 6 and Data Set 4 in Table 5 N"6
h E(B ) E(B ) E(B ) E(B ) E(B ) E(B ) E(B )
N"9
N"12
Algo
Sim
% Error
Algo
Sim
% Error
Algo
Sim
% Error
0.698 1.744 4.094 1.906 2.881 2.881 1.379 1.379
0.710 1.888 4.121 1.879 2.420 2.336 1.691 1.778
! 1.69 !7.63 !0.66 1.44 19.05 23.33 !18.45 !22.44
0.778 2.625 6.044 2.955 4.289 4.289 2.091 2.091
0.789 2.825 6.170 2.830 3.642 3.461 2.533 2.714
!1.39 !7.08 !2.04 4.42 17.76 23.92 !17.45 !22.96
0.820 3.566 7.898 4.101 5.681 5.681 2.750 2.750
0.834 3.788 8.178 3.822 4.810 4.623 3.402 3.589
!1.68 !5.86 !3.42 7.30 18.11 22.89 !19.17 !23.38
Fig. 4. A six-machine network model. (Machines 4}6 are leaf nodes).
additional state variables to capture the additional information. As the partition gets "ner, the model becomes more accurate. 6.2. Adding tunable parameters Estimates (4.1)}(4.3) of the aggregate transition rates can be augmented with additional parameters in order to meet additional conditions on the system. This is illustrated by one approach to ensure that Eqs. (2.5) and (2.7) hold: Expression (4.2) for the estimated Kanban arrival rate cD* as seen by machine f, for f3INT, can be modi"ed by incorporating an additional parameter a . This is justi"ed by the recognition that D
1036
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
Fig. 5. A "ve-machine network model. (Machines 4 and 5 are leaf nodes).
Fig. 6. A six-machine network model. (Machines 3}5 are leaf nodes).
the approximation Eq. (3.5) is not exact. The new parameter helps determine cD* and hence h , and D it should be chosen so that h "h . Satisfaction of this condition for all fOr is Eq. (2.5). D P Similarly, expression (4.1) for the e!ective service rate kD* at machine f3INT can be modi"ed by employing an additional parameter b . This is justi"ed by the recognition that approximation (3.1) D , and it should be chosen so is not exact. This new parameter helps determine kD* and hence En D1D that En "Em !Em . D1D 1D D Satisfaction of this condition for all f3 INT will ensure that the path-independence condition (2.7) holds.
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1037
6.3. Partitions with multiple machines Simulation shows the existence of a `#oating bottlenecka, where completion of an assembly at a slow machine triggers a quick succession of other assembly operations downstream, with this cascade stopping when the platoon of parts reaches a new downstream bottleneck. The aggregate model of f fails to capture this dynamic correlated behavior of several machines, which may be the source of the large errors in computing mean queue lengths, a common failing of `single-nodea decompositions. The appropriate remedy is a partition that simultaneously encompasses several assembly machines. This, unfortunately, leads to an explosion of the number of the states. A coarser partition of the Cartesian-product state space seems appropriate. Computational results for blocking networks [13] suggest that there are diminishing marginal bene"ts from incorporating more than a few machines in the aggregate model. 6.4. Hybrid schemes Hazra [1,7] describes a di!erent aggregation approach for evaluating assembly networks, and surveys other available techniques. Each technique has its own strengths and weaknesses. An open issue is how to construct hybrid approaches that combine the bene"ts of several pre-existing schemes. For example, Hazra's [1] scheme has only 3% error in computing mean queue lengths, while the present A/D scheme has more accurate throughputs. Could the aggregate service rates from one scheme be used as inputs to the other? The goal would be to achieve the advantages of both: only 1% error in throughput and only 3% error in mean queue length. 6.5. Triangular decomposition of the xxed-point problem There is structure in the "xed-point equations that has not yet been exploited. The in#uence diagram for the network in Fig. 1 is given by Fig. 7. Fig. 7 shows the interactions in the "xed-point system, and mimics the physical topology. At machine f, the service rate information kD is passed downstream to the successor while the Kanban arrival rate information cD is passed upstream to all the predecessors. For example, the model at machine 2 receives k information from the machine-4 model, k information from the machine-5
Fig. 7. In#uence diagram for the assembly network in Fig. 1. Note: for f3¸"+7, 8, 9, 10, 11, 12,, kD are already known exactly while cD are not needed.
1038
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
model, and c information from the machine-1 model. The model at machine 2 then generates c information, which it sends to the machine-4 model, c information, which it sends to the machine-5 model, and k information, which it sends to the machine-1 model. Each machine model interacts with all its immediate neighbors, and with no one else, and these information #ows are always bidirectional. The resulting equations are block triangular. For example, the Kanban arrival information c passed to machine 2 contains all the information needed to compute the behavior at the subtree consisting of machines 2 and 4}10. In particular k and k will be implicit functions of c, so that the entire subtree can be collapsed into a single ,ctitious machine 2I. This provides the justi"cation for an aggregation approach whereby the subtrees are merged sequentially, similar to the merging in [7]. The critical feature in any such merging scheme, of course, is making accurate estimates of the queue-dependent service rate kI of the "ctitious machine 2I. Hazra's results show that this approach is highly promising.
7. Summary The performance evaluation of a closed Kanban-controlled assembly system is accomplished by a new heuristic version of the exact aggregation}disaggregation theory for "nite continuous-time Markov chains. It provides accurate estimates of the throughput. Because the approximation is theory based, it provides a framework for further model development, with some possible extensions described in the paper. The approximation has the novel feature of doing simultaneous multiple partitions of the state space, in such a way that the associated aggregate transition rates are mutually consistent. The methodology is a novel approach towards extending aggregation ideas to fork-and-join queueing networks. It provides fast and accurate throughput estimates for complex production systems, and insights into the behavior of these systems. We use the results of this methodology to show how the throughput of an assembly system varies with changes in the parts #ows, in the total number of Kanbans, or in the processing rate at the individual stations. The speed at which we obtain these results make it a practical tool for the analysis of key design trade-o!s for complex manufacturing systems. Anticipated use of the model includes balancing of routes, speeds and topologies.
Acknowledgements Parts of this paper were completed while P. Schweitzer was a Lady Davis Fellow at Hebrew University in Jerusalem, Israel. He thanks the University and Lady Davis Trust for many courtesies extended.
References [1] Hazra J, Seidmann A. Performance evaluation of closed tree-structured assembly systems. IIE Transactions 1996; 28:591}9.
J. Hazra et al. / Computers & Operations Research 26 (1999) 1015}1039
1039
[2] Takahashi Y. A lumping method for numerical calculations of stationary distributions of Markov chains. Research Report No. B-18, Dept. of Information Sciences, Tokyo Institute of Technology, Tokyo, Japan, 1975. [3] Schweitzer PJ. Aggregation methods for large Markov chains. In: Iazeolla G, Courtois PJ, Hordijk A. editors. Mathematical computer performance and reliability. Amsterdam: North-Holland, 1984:275}86. [4] Schweitzer PJ. A Survey of aggregation}disaggregation in large Markov chains. In: Stewart WJ, editor. Numerical solution of Markov chains, New York: Marcel Dekker, 1991:63}88. [5] Dallery Y, Liu Z, Towsley D, Equivalence, reversibility, Symmetry and concavity properties in fork-join queueing networks with blocking. Journal of Association for Computing Machinery 1994;41:903}42. [6] Viswanadham N, Narahari Y. Performance modeling of automated manufacturing systems. Englewood Cli!s, NJ; Prentice-Hall, 1992. Hoop MJ, Spearman ML. Factory physics, foundation of manufacturing management. Homewood, IL: Irwin, 1996. [7] Hazra J. Performance evaluation of Kanban-controlled assembly networks. PhD Thesis, William E. Simon Graduate School of Business Administration, Rochester, NY: University of Rochester, 1995. [8] Gershwin SB, Manufacturing systems engineering. Englewood Cli!s, NJ: Prentice-Hall, 1994, Buzacott JA, Shanthikumar JG. Stochastic models of manufacturing systems. S. Englewood Cli!s, NJ: Nahmias, Prentice-Hall, 1993 Production and operations analysis. Homewood, IL: Irwin, 1997. [9] Chandrasekhar Rao P, Suri R. Performance analysis of an assembly station with input from multiple fabrication lines. Working paper Dec. 1996, Department of Industrial Engineering, University of Wisconsin-Madison, Madison, WI. [10] Chandrasekhar Rao P, Suri R. Approximate queueing network models for closed fabrication!assembly systems } Part I: Single level systems. Production and Operations Management 1998;3:244}75. [11] Baynat B, Bruno, Dallery Y. Approximate techniques for general closed queueing networks with subnetworks. EJOR 1992; 69(2): 250}64. [12] Schweitzer PJ. Aggregative modelling of queueing networks. In: Bonatti M, editor. Teletra$c science for new cost-e!ective systems, networks and services. Amsterdam: Elsevier, North-Holland, 1989:1522}8. [13] Takahashi Y. Aggregate approximation for acyclic queueing networks with communication blocking. In: Perros HG, Altiok T, editors. Queueing networks with blocking. Amsterdam: Elsevier, 1989:33}46. [14] Dimascolo M, Frein Y, Dallery Y. An analytical method for performance evaluation of Kanban controlled production systems. Operations Research 1996;44(1):50}65.
Jishnu Hazra received his B.Tech (Mech. Engg.) from IIT, Delhi and M.S. and Ph.D. from the University of Rochester in Operations Management. His research interests are in Supply Chain Management and JIT Manufacturing. He has conducted Executive Development Programs on Supply Chain Management for TATA}IBM executives and other companies. He is also involved in consulting and research projects in the areas of R&D Project Planning and Logistics. Paul J. Schweitzer is Professor of Business Administration at the William E. Simon Graduate School of Business Administration of the University of Rochester. His interests include military operations research, performance models of manufacturing and computer systems, Markovian decision processes (MDPs), management science and information technology. He has done seminal work on value-iteration for MDPs, perturbation of Markov chains, mean value analysis of queueing networks including #exible manufacturing systems, and aggregation techniques. He is the author of over 100 research papers and reports. Abraham Seidmann is the Xerox Professor and Areas Coordinator of Computers and Information Systems, Management Science and Operations Management at the William E. Simon Graduate School of Business Administration at the University of Rochester. Professor Seidmann is the author of over 75 research articles which appear in many of the leading scienti"c journals and is a Department Editor on Interdisciplinary Management Research and Applications in Management Science. His current research and consulting activities include Information Systems, Business Process Reengineering, Project Management and Optimal Resource Allocation, Strategic Manufacturing Systems, Information Economics, Healthcare Management, Stochastic Processes and Performance Modeling for Capacity Planning of Assembly Systems.