Analysis of assembly systems controlled with kanbans

European Journal of Operational Research 166 (2005) 310–336 www.elsevier.com/locate/dsw

Production, Manufacturing and Logistics

Analysis of assembly systems controlled with kanbans Andrea Matta b

a,*

, Yves Dallery

b,1

, Maria Di Mascolo

c,2

a Dipartimento di Meccanica, Politecnico di Milano, via Bonardi 9, 20133 Milan, Italy Laboratoire Productique-Logistique, Ecole Centrale Paris, Grande Voie des Vignes 92295, Chatenay-Malabry Cedex, France c Laboratoire d’Automatique de Grenoble, 38402 Saint-Martin d’H
eres, Grenoble, France

Received 27 May 2002; accepted 4 March 2004 Available online 11 May 2004

Abstract This paper deals with the problem of evaluating the performance of assembly systems managed with kanbans. In particular we analyze kanban systems functioning with different control policies depending on how kanbans are released in assembly stages. Simultaneous and independent releases of kanbans are considered in the analysis as alternative control policies to be used in assembly systems. Queuing network techniques are used to calculate the major system performance measures such as throughput of the system, percentage of satisfied orders, average level of finished products and average time delay of unsatisfied orders. Simultaneous and independent kanban assembly systems are compared in some numerical cases. The differences, in terms of performance, between the two alternative kanban control policies are evaluated so that the proper control policy can be selected in the design phase of the system. Simulation results are also reported to validate the approximation of the analytical method used to calculate the performance of assembly kanban systems. Ó 2004 Elsevier B.V. All rights reserved. Keywords: Queuing; Control; Kanban; Production systems; Assembly lines

1. Introduction Kanban controlled production systems have received much attention in the past because of the ease in which they can implement the pull control policy using a make to stock production (see [7,11,16,17] for details and references). The production system is decomposed into several stages, each of which corresponds to a subpart of the system, where production is controlled by way of cards (called kanbans), acting as production orders. The design of the control system consists in determining the number of kanbans that must be associated with each stage. In order to achieve this goal, a performance evaluation of the system is

*

Corresponding author. Tel.: +39-2-23994925; fax: +39-2-70638377. E-mail addresses: [email protected] (A. Matta), [email protected] (Y. Dallery), [email protected] (M. Di Mascolo). 1 Tel.: +33-1-41131533; fax: +33-1-41131272. 2 Tel.: +33-4-76826416; fax: +33-4-76826388.

0377-2217/$ - see front matter Ó 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2003.09.035

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

311

required, for given values of the number of kanbans: this is then used as a building block in an optimization procedure. The paper focuses on assembly kanban systems functioning with different control policies, depending on whether the releases of kanbans at assembly stages take place independently or simultaneously. Indeed, system performance depends on the type of control policy that is adopted: up to now there is no study that points out quantitatively the main differences. The kanban control system with simultaneous release was considered in [9,10], and the kanban control system with independent release was introduced in [6]. The paper has two goals. First a new approximate analytical method is proposed for evaluating the performance of assembly systems with independent release of kanbans. This method is mainly based on product-form approximation [1,2] and multi-class aggregation techniques [3] developed by Baynat and Dallery. Our purpose here is to extend the multi-class approximation technique [3,4] to analyze the performance of assembly systems. The second goal of the paper is to compare the two different control policies for releasing kanbans in assembly systems: a numerical comparison has been carried out for identifying the most suitable application areas of the analyzed control policies. The paper is structured as follows. Section 2 introduces production systems managed with kanbans, Section 3 contains a description of assembly kanban systems. In Section 4 the methodology used to evaluate the system performance is briefly explained and numerical results are reported in Section 5. Finally, conclusions are drawn in Section 6. Technical details of the proposed analytical method are also described in Appendix A of the paper.

2. Kanban mechanism A simple production system is shown in Fig. 1 where ellipses and triangles represent the production stages and the buffers of the system respectively. We assume that parts spend a certain amount of time in each production stage in order to complete their process cycle; this time is not deterministic and mainly depends on the specific operation that is performed on parts (e.g. casting, cheap removal, packaging, etc.). Parts follow a linear path starting from the first stage until the last stage is reached; then finished products wait in the last buffer (i.e. the buffer N in Fig. 1) until a request for one part is sent by customers. Each stage can contain one or more machines; for instance, a stage can be a job shop made up of turning machines, drillings, grindings, etc. in which several cheap removal operations are performed on the parts. Buffers separate production phases in order to decrease interdependencies between the different stages. The part flow entering each stage is controlled so that the amount of work in process (WIP) in the system is limited and the related inventory costs are reduced. In particular, parts have to receive an authorization order before entering process stages. The arrival of an authorization to access a stage of the system is modeled by the arrival of a kanban, normally a card in real shop floors. Thus, the number of kanbans, which represent the authorization orders, controls the flow of parts in the system. Let us denote by FPi the buffer containing products that have been manufactured by stage i (where i ¼ 1; . . . ; N ) and by ki the number of kanbans at stage i (Fig. 2 shows the links between a stage i of a kanban system with stages i
1 and i þ 1, using queuing network modeling). Parts in buffer FPi wait for a production order to be released by the system controller, represented by entities in queue Diþ1 . Since i-stage parts enter stage i, each part has attached an i-stage kanban (i.e. the part

Raw Parts

Buffer 0

Stage 1

Buffer 1

…

Buffer N-1

Stage N

Fig. 1. Production system with N stages.

Buffer N

Finished products

312

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336 FPi-1 (i-1)-stage products, (i-1)-stage kanbans

FPi

Si i-stage products, (i+1)-stage kanbans

Di

Di+1 ki

ki-1

k i+1

Fig. 2. Stage i of a kanban system.

has the authorization to flow in stage i). Kanbans in buffer Di , which represent the production orders of stage i, wait for the arrival of ði
1Þ-stage finished parts while kanbans in buffer Diþ1 , which represent the production orders of stage i þ 1, wait for the arrival of i-stage finished parts. We now describe the flow of physical parts in stage i of the system. To proceed to stage i, each (i
1)stage finished product has to find at least one unit of kanban in buffer Di . If a part and a kanban are contained in FPi
1 and Di respectively, the part proceeds together with the kanban into stage i where it will continue its process cycle; we say the part synchronizes with the kanban before entering stage i. After synchronization occurs, the (i
1)-stage kanban is detached from the physical part and is carried over to the previous stage. If buffer Di is empty, the (i
1)-stage finished part waits in FPi
1 until a new production order is released; this corresponds to the case in which many orders have already been released and are still being processed in stage i. A general serial system with N stages controlled with kanbans is modeled with the queuing network shown in Fig. 3. The Kanban control mechanism provides close coordination between stages. The information that the number of customer orders has increased is not immediately transmitted from the last stage to upstream stages. This information delay can reduce system WIP as well as its corresponding inventory costs. The reaction time and, therefore, the delay of transmission of consumption of demand to the upstream stages of the system depend on the number of kanbans at each stage: the larger the number of kanbans, the quicker the system will react. Refer to [5,13,16,17] for a comparison of the kanban mechanism with other pull management policies such as base stock, generalized kanban [12,14] and extended kanban. Among the methods proposed in the literature for analytical analysis of serial kanban systems (see [18,19]), that presented in [11] appears to be of special interest since it can handle stages consisting of any number of machines and is fairly accurate.

3. Assembly kanban systems In assembly systems, different part flows join at some points along the flow path; these joining points are assembling stages where two or more components are combined to form a whole product. Let us consider

Raw Parts

FP0

FP1

FP2

S1

D1

D2 k1

FPN-1

S2

…

DN

D3 k2

FPN

SN

k3

kN-1

DN+1 kN

Fig. 3. N -stage production system controlled with kanbans.

Demand

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

313

an assembly system with N stages where Ri
1 is the number of types of components leaving stage i
1 and assembled at stage i. System behavior is similar to that of a kanban system described in the previous section except for the assembling stations in which components are combined. The queuing network shown in Fig. 4 represents an assembly system with N stages and two components. We denote by ki;r the number of kanbans at stage i related to the component of type r (where r ¼ 1; . . . ; Ri ). In the example shown in Fig. 4, the parameter Ri
1 is equal to 2; raw materials enter the first stage of the system both for components of types 1 and 2. Components are processed until they leave stage i where one unit of type 1 is assembled with a unit of type 2; then, the assembled part continues its processing cycle until the final operations are completed at stage SN . However, before entering stage i, components have to wait for authorization. Let us enter into more detail. Each unit of an (i
1)-stage component of type r, which has completed stage i
1, waits in queues FPi
1;r until the production order of stage i is available, i.e. at least one i-stage kanban is in queue Di . An assembly operation can be carried out only if both queues FPi
1;r and Di contain at least one unit. We model the joining operation with a synchronization station in which Ri
1 different components are combined at stage i. A synchronization station between stages Si
1 and Si is denoted by Ji
1;i . The way in which kanbans from the previous stage are released in the assembling station depends on the specific rule adopted by the system. A first simple management rule allows release of kanbans from the previous stage only after the assembly operation can start (see Figs. 4 and 5). With such a rule the demand information is transmitted at the same time to the Ri
1 upstream stages that produce the components, i.e. release of kanbans is simultaneous. Even if a component leaves stage i
1 before the others, its corresponding kanban is released only when the assembly operation can start. Therefore, the transmission of demand information to upstream stages can be delayed in the case of some components, in particular for the first ones to enter the station, while the release of kanbans is not delayed only for the last part reaching the synchronization station. These types of systems are known as simultaneous kanban systems [10]. Performance of simultaneous kanban systems can be evaluated by using simulation or by applying the method proposed by Di Mascolo and Dallery in [10] and based on previous works [2–4,8,11]. A second management rule does not allow any delay in releasing the kanbans. Kanbans of the previous stage are released immediately after the arrival of components and the authorization to access the next stage. In particular, one unit of r-type kanbans of (i
1)-stage is released immediately if at least one

ki-2,1

k1,1 D1,1

ki-1,1

D i-2,1 FPi-1,1 Si-1,1

Raw parts

… FPi-2,1

FP0,1

Si Raw parts

… FP0,2

FPi FPi-2,2

k1,2

Di-1,2

FPi

FPN-1

Di+1

DN

FP N

Si-1,2 FPi-1,2

D1,2

Finished Products

SN

…

ki-1,2

ki+1

kN-1

ki-2,2 Di ki

Fig. 4. N -stage assembly system with two components.

D N+1 kN Demand

314

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336 FPi-1,1

Component 1

FPi-1,2

Component 2

...

... FPi-1,R

Assembled part and i-stage kanban

Component R Di k i-1,1 … k i-1,R

Demand ki

Fig. 5. Simultaneous assembly synchronization station.

component is in queue FPi
1;r and one production order is in queue Di . Dallery and Di Mascolo call these types of systems independent kanban systems [10]. When a component leaves stage i, its kanban is released also if the assembly operation cannot start; the only condition for releasing the kanban is that at least one rtype component is in queue FPi
1;r and one production order is in queue Di . Therefore, transmission of information on demand consumption to the upstream stage is not delayed. We model the behavior of independent assembly synchronization stations with a set of synchronization stations (see Fig. 6). The behavior of the station is as follows. When a new authorization order reaches the synchronization station from stage i þ 1, instances of kanbans are generated so that the same information is sent to different Ri queues. In the example shown in Fig. 6, one kanban enters the station and generates one instance so that two kanbans enter queues D1 and D2 respectively. Queues FP1 and FP2 contain the i-stage components that reach the station; each of these components has an attached kanban of stage i. Kanbans attached to 1st-type components can leave the station only if at least one authorization code is in queue D1 and one unit is in queue FP1 . If this condition is satisfied, the i-stage kanban is detached from the component and is sent to the previous stage; in addition an (i þ 1)-stage kanban is attached to the 1st-type component and enters the queue C1 where it waits to synchronize with a unit of queue C2 . In the same way, kanbans attached to 2nd-type components can leave the station only if at least one authorization code is in queue D2 and one unit is in queue FP2 . If this condition is satisfied, the i-stage kanban of type 2 is detached

ki,1

FP1 Component 1

C1 D1

Demand

D2

Assembled component

C2 and (i+1)-stage kanban Component 2

FP2 k i,2

Fig. 6. Independent assembly synchronization station.

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

315

from the component and is sent to the previous stage. Also in this case an (i þ 1)-stage kanban is attached to the 2nd-type component and enters the queue C2 where it waits to synchronize with a unit of queue C1 . Synchronization between C1 and C2 completes the assembly operation. Performance of independent kanban systems can be evaluated by means of simulation or using the new method proposed in this paper. In the next section we briefly explain the analytical approximation technique used in this paper for evaluating the performance of the described assembly systems. In particular the proposed methodology is similar to that developed in [10] for simultaneous kanban systems except for analysis of assembly synchronization station in an isolation mode, for analysis of which a new method based on the multi-class aggregation technique [3] has been proposed in this paper. Then numerical analysis of different types of kanban control rules in assembly systems is reported in Section 5.

4. Performance evaluation model 4.1. Outline of the method Each server of the queuing network representing the assembly system can model the behavior of production stages or that of synchronization stations. If production stages contain more than one machine, the model can be further enlarged by adding nodes to the network so that the behavior of the stages is properly represented. Synchronization station Ji;iþ1 combines Ri finished components of stage i and a free kanban of stage i þ 1. In particular (i þ 1)-stage kanbans enter station Ji;iþ1 and synchronize with the Ri components waiting to enter stage i þ 1; if at least one buffer FPi;r is empty, kanbans in Diþ1 wait in the queue. The arrival of physical parts or components in synchronization stations can be viewed as a waiting time for production orders or kanbans. For instance, kanbans of stage i þ 1 wait in station Ji;iþ1 until all the parts necessary for the synchronization have arrived. As a consequence, it is possible to model synchronization stations as single servers where kanbans arrive and wait in the queue for components to arrive. The entire system is modeled as a multi-class queuing network L where kanbans are the populations of the network, stages (or machines and buffers) and synchronization stations are the servers. The total sum of servers of the network includes both real machines and synchronization stations. Note that the behavior of the different network classes is not independent because they interact at synchronization stations of different adjacent stages. Indeed, when a kanban synchronizes with parts, the kanbans of previous stages are detached from the parts and released to the upstream phase. Therefore, synchronization of kanbans of stage i þ 1 at station Ji;iþ1 affects the release of kanbans of stage i; again synchronization of kanbans of stage i at station Ji
1;i affects the release of kanbans of type i
1 and so on. In order to evaluate the performance of network L, we use a general methodology based on a productform approximation technique [1,2]. The idea of the method is to decompose the multi-class system L into a set of single-class closed queuing networks Ni;r , each one related to a specific production stage and type of kanban. The method associates N single class product-form networks, each one modeling the behavior of a particular network Ni;r , to the original multi-class queuing system. Each sub-network is visited by only one class of customers (we denote by Sr the set of servers visited by class r). In particular, kanbans of type i; r, i.e. related to stage i and component r, can visit machines of sub-network Ni;r . The total number of single-class queuing networks is equal to the number of different types of kanbans in the system. Sub-networks are closed, i.e. the number of kanbans flowing in sub-network Ni;r is constant and equal to ki;r . Since each single-class is a Gordon–Newell network, the stationary solution of the network Ni;r is known and has the following form [15]:

316

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

P ðnÞ ¼

1

" k i;r Y Y

Gðmi;r ; ki;r Þ m2Sr

# Vm;r ; l ðn Þ nm;r ¼1 m;r m;r

ð1Þ

where n is a vector representing the global state of the system, i.e. the number of customers at each server, P ðnÞ is the steady state probability of the analyzed network, mi;r is the total number of stations in the subnetwork Ni;r (synchronization stations are included), Sr is the set of servers visited by class r, and G is a normalization constant. Parameters Vm;r and lm;r ðnm;r Þ are respectively the visit ratios and the service rates of network servers. Thus, the performance of the multi-class in the original network is approximated by the performance of single class networks. The methodology is obviously an approximation since it assumes independence among the different classes. However, each single class network has stations with load-dependent service rates. This allows consideration in a sub-network, in particular in evaluation of load dependent service rates, of the effect of the rest of the system. For this purpose, the service rates lm;r ðnm;r Þ of load dependent servers have to be properly calculated. In particular, service rates are evaluated by taking into account the analysis of servers in an isolation mode. The analysis of servers in an isolation mode, i.e. the stations functioning independently from the rest of the system, provides the throughputs of the stations considering more than one class so that interactions among them are partially captured. Ideally, the service rates lm;r ðnm;r Þ of network servers should be equal to the conditional throughputs mm;r ðnm;r Þ of class r customers at the corresponding station m in the original network. Since throughputs of servers in isolation represent a good approximation of conditional throughputs, we analyze the servers in isolation to approximate the service rates of load dependent servers. The method is iterative. The first step of the methodology is to calculate the average arrival of customers entering each node of the network. This task is fairly straightforward because single-class closed networks have to be solved using traditional techniques. Once the average arrival rates of single servers are known, servers are analyzed in an isolation mode so that stationary state probabilities of the subsystem can be calculated as well as station throughput. Throughputs of single servers analyzed in isolation mode are assumed (note that an approximation is introduced in this step) to be the service rates of load dependent servers of sub-networks Ni;r . The method is iterative until the unknown parameters lm;r ðnm;r Þ converge to the stationary values; in particular service rates are solutions of a fixed-point problem and the number of iterations to achieve convergence is usually very reasonable. Refer to [1,2] for more details on the algorithm. After the proper values of load dependent servers have been calculated, performance indicators can be easily evaluated by solving the single-class closed networks. 4.2. Analysis of assembly kanban synchronization stations in isolation In the overall algorithm, a performance evaluation tool must be used for calculation of throughput of servers in isolation. We describe in this section the analytical methods used to analyze simultaneous and independent assembly synchronization stations. For classical queuing servers (e.g. M/M/1, M/M/m, M/Er/ 1, etc.) solutions are already well known in technical literature and are not dealt with in this paper. 4.2.1. Simultaneous assembly synchronization stations A model of a simultaneous assembly synchronization station is shown in Fig. 7. The synchronization station Ji;iþ1 is made up of Ri þ 1 different queues crossed by Ri þ 1 different classes corresponding to Ri components and one production order. The first Ri queues represent buffers FPi;r containing components that have to be assembled while the last queue contains production orders of stage i þ 1. One unit of each class synchronizes simultaneously with one unit of the demand class. We denote by ki;r ðni;r Þ the average arrival rate, calculated by solving the closed single-class networks, of components into queue FPi;r and by

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

Component 1

Component 2

317

FPi,1

λi,1(ni,1)

FPi,2

λ i,2 (n i,2)

...

... Component R

λi,R (ni,R)

λi,D (ni,d) k i,1 … ki,R

FPi,R

Assembled part and i-stage kanban

D i+1

Demand k i+1

Fig. 7. Simultaneous assembly synchronization station Ji;iþ1 .

ni;r the number of parts that are contained in that queue. Note that the arrival of parts in the synchronization station depends on the number of parts in the corresponding queues. This is due to the fact that the number of kanbans of each stage is constant and limits the number of parts circulating in that phase of the system. Furthermore, also the last queue, denoted by Diþ1 has an average input kDiþ1 ðnDiþ1 Þ of kanbans that depends on the buffer level nDiþ1 . Thus, the system is a multi-class synchronization station with state dependent arrival rates. Synchronization can take place only if at least one unit is present in each queue. After synchronization is complete, the assembled unit is sent to the next stage while kanbans attached to the components are simultaneously sent to the previous stage. Synchronization time is assumed to be deterministic and equal to zero. The system behavior can be modeled as a Markovian process where the vector n ¼ ðni;1 ; ni;2 ; . . . ; ni;R ; nDiþ1 Þ represents all its possible states. In order to evaluate the throughput of the synchronization station, it is necessary to calculate the stationary probabilities of each state of the system. An exact solution can be found by applying traditional numerical methods. However, this approach is feasible only for very small systems (after 5 queues with 20 customers for each queue, the computational time is quite long). For large systems, the total number of states of the Markov chain becomes too big and the approach may be unfeasible for practical applications. Approximate solutions are available in the literature. In particular, Baynat and Dallery [3] proposed an approximate method to calculate the performance of multi-class synchronization stations with state-dependent arrival rates. The technique consists in representing the simultaneous synchronization station in a set of smaller synchronization stations or blocks easier to analyze than the original system (see Fig. 8). Each small synchronization station models the behavior of the entire synchronization station and is made up of two queues: the first queue represents the behavior of class r in queue FPi;r while the second queue represents aggregately that of the other classes. In a decomposition of a synchronization station with Ri þ 1 classes, the number of blocks is equal to Ri þ 1. We use the notation 2qSS to refer to these blocks. Parameters of blocks (i.e. arrival rates and queue capacities) must be appropriately defined so that the throughput of parts leaving the station is equal for all classes. A set of equations concerning the equality of throughputs among the blocks and the number of parts contained in the system allows evaluation of the 2qSS parameters. Thus, there is a total of Ri þ 1 two-queue synchronization stations that globally model the entire synchronization station. The complexity of the problem is therefore reduced since an exact closed solution of a 2qSS is already available in the literature [8]. The method is implemented in an iterative algorithm solving a fixed point problem. Starting from an initial solution, the algorithm changes the arrival rates of each 2qSS

318

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336 FPi, 1

Class 1 Aggregated Class FPi,1

Class 1 Class 2

... Class R

FPi,,1

Class 2 FPi,2

Synchronized part

Aggregated Class

...

...

FPi,,1

FPi,R Synchronized

part Di+1

ki,1 … ki,R

Synchronized part

Demand

Class R Synchronized part

Aggregated Class Di+1

ki+1

Class R+1 Aggregated Class

Synchronized part

Fig. 8. Multi-class aggregation [3].

until their throughputs converge to the unknown value. In particular a 2qSS is solved for each type of class that visit the system (in total Ri þ 1) on the basis of reasonable starting assumptions on the unknown arrival rates. Throughputs of each 2qSS are then calculated and the values of arrival rates for each 2qSS are modified so that their throughput is equal to the average of the throughputs calculated in the previous iteration. The algorithm continues until throughputs of the Ri þ 1 blocks converge to the same value. The advantage of the technique is its lack of complexity, since it requires Ri þ 1 solutions of 2qSS instead of solving one synchronization station with Ri þ 1 queues. In this paper, Baynat and Dallery’s method [3], successively adopted in [4,10], is used to calculate the performance of simultaneous kanban systems. 4.2.2. Independent assembly synchronization stations Let us now analyze independent assembly synchronization stations. An overview of the new analytical method for evaluating the performance of independent kanban systems is provided in this subsection while the technical details are reported in Appendix A. For simplicity’s sake, we describe the case in which two types of components have to be synchronized with kanbans of the next stage, but the modeling is general enough to be extended to the case of Ri components. Queues FPi;1 and FPi;2 contain components of types 1 and 2 respectively that have been produced up to stage i. Each component is assigned a kanban of stage i; queues Diþ1;1 and Diþ1;2 contain the authorization codes, i.e. kanbans, to allow components to enter the next stage i þ 1. We denote by ki;r ðni;r Þ (where r ¼ 1; 2) the average arrival rate of components into queue FPi;r and by ni;r the number of parts that are contained in that queue. Kanbans of type r are released if at least one unit is contained in queues FPi;r and Diþ1;r ; that is, the release of kanbans of type r (with r ¼ 1; . . . ; Ri ) is independent and does not depend on the release of the other types of kanbans. The described independent assembly synchronization station is represented by the network shown in Fig. 9. Kanbans derived from stage i þ 1 are duplicated so that one unit goes into queue D1 and the duplicated one into queue D2 (index i is omitted for simplicity). Components enter the system and wait for the production order in queues FP1 (for components of type 1) or FP2 (for components of type 2); index i is

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

319

ki,1

FP1

λ1(nFP1)

C1

λd(nD)

D1 D2

Assembled component

C2 and (i+1)-stage kanban λ2(nFP2) FP2 ki,2

Fig. 9. Independent assembly synchronization station Ji;iþ1 .

omitted for simplicity. Components from stage i synchronize with those contained in queues D1 or D2 . Previous stage kanbans leave the system after crossing the synchronization stations made up of queues FP and D. Assembled parts leave the system only after the last synchronization station of queues C1 and C2 (index i is omitted for simplicity). The state of the system can be described by the vector n ¼ ðnFP1 ; nFP2 ; nD1 ; nD2 ; nC1 ; nC2 Þ representing the levels of queues FP1 , FP2 , D1 , D2 , C1 and C2 . These queues are not independent. The sum of buffer levels in queues D1 and C1 must be equal to the number of kanbans contained in the original queue Diþ1 . Likewise, the sum of buffer levels in queues D2 and C2 is equal to the number of kanbans contained in the original queue Diþ1 Therefore, we can write the following quantities: nD1 þ nC1 ¼ nD2 þ nC2 ¼ nDiþ1 : Also, the total number of components of type 1 in the station is equal to the sum of buffer levels of queues FP1 and C1 . Alternatively, the total number of components of type 2 in the station is equal to the sum of buffer levels of queues FP2 and C2 . Thus, the described system is equivalent to a single synchronization station made up of only three equivalent queues (see Fig. 10). Let us denote the queues of the equivalent system by the notations E1 ; E2 ; ED , their corresponding levels by ne1 ; ne2 , and neD , and their arrival rates by ke1 ðne1 ; neD Þ, ke2 ðne2 ; neD Þ and keD ðneD Þ (for simplicity, we omit index i standing for the number of the stage). In order to calculate the performance of station Ji;iþ1 it is sufficient to analyze only one synchronization station instead of the whole system. The vector ne ¼ ðne1 ; ne2 ; neD Þ describes the states of the equivalent system. If R different types of components enter station Ji;iþ1 the equivalent vector describing the system becomes  ne ¼ ðne1 ; . . . ; ner ; . . . ; neR ; neD Þ. To move from the original system to the equivalent one, the following transformations can be used:

λ1e (n1e , n De ) e

e

λ2 (n 2 , nD )

e

e

λD (n De)

E1

E2

ED

Assembled part and (i+1)-stage kanban

Fig. 10. Equivalent model of the independent assembly synchronization station Ji;iþ1 .

320

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

8 e n1 ¼ nFP1 þ nC1 > > > e > > n2 ¼ nFP2 þ nC2 > > > < ne ¼ n þ n ¼ n þ n D1 C1 D2 C D
2 e e  e e e e > k1 ðn1 ; nD Þ ¼ k1 ðn1
min n1 ; nD Þ > >
e e > e e > e e > k > 2 ðn2 ; nD Þ ¼ k2 ðn2
min n2 ; nD Þ > : e e kD ðnD Þ ¼ kD ðneD Þ;

with with with with with

ne1 ¼ 0; 1; . . . ; min fki;1 ; kiþ1 g þ ki;1 ¼ N1e ; ne2 ¼ 0; 1; . . . ; min fki;2 ; kiþ1 g þ ki;2 ¼ N2e ; neD ¼ 0; 1; . . . ; kiþ1 ¼ NDe ;
 max 0; ne1
ki;1 6 neD 6 kiþ1 ;
e  max 0; n2
ki;2 6 neD 6 kiþ1 ;

ð2Þ

where ki;1 and ki;2 represent the maximum number of components of types 1 and 2 respectively that can enter system Ji;iþ1 ; kiþ1 is the maximum number of arrivals of production orders. The average arrival rates of entities in the equivalent queues are state dependent. In particular, the arrival rate into queue E1 depends on the level ne1 of its queue but also on the level neD of queue ED . Similarly, for components of type 2, the arrival rate into queue E2 depends on the level ne2 of its queue and on the level neD of queue ED . Finally, the arrival rate into queue ED depends only on the level neD of its queue. In order to evaluate the throughput of the station, it is necessary to analyze the Markov chain with state dependent transition probabilities of the equivalent system. Also in this case, an exact solution can be found by applying traditional numerical methods. However, this approach is feasible only for very small systems. For large systems, approximate techniques must be used to evaluate the system performance. In this case it is not possible to apply the multi-class aggregation technique used by Baynat and Dallery [3]. The reason is that the arrival rates ker ðner ; neD Þ in the equivalent systems depend on queue levels ner and neD . Indeed, if Baynat and Dallery’s technique was used, it would not be possible to aggregate in an equivalent class the behavior of the class related to the demand. This is because of the statedependence of arrival rates. Queues Er (if class r is analyzed) and ED have to be analyzed explicitly so that arrival rates ker ðner ; neD Þ can be calculated. However, the 2qSS blocks include the behavior of the class related to the demand in an aggregated class. It is not possible in this way to know the state of queue ED since its behavior is modeled in the aggregated queue class; as a consequence arrival rates cannot be calculated using this technique. In this paper we develop an extension of Baynat and Dallery’s technique [3]. This extension models the behavior of the synchronization station with a set of small synchronization stations or blocks. Blocks are made up of three queues: the first queue represents the behavior of class r in queue FPi;r , the second queue represents that of class R þ 1 (i.e. the class related with the demand), while the third queue represents aggregately that of the other classes. We use the notation 3qSS to refer to these blocks. Similarly with the approach adopted by Baynat and Dallery, parameters of blocks (i.e. arrival rates and queues capacities) must be selected appropriately so that the throughput of parts leaving the station is equal for all the classes. A set of equations concerning the equality of throughputs among the blocks and the number of parts contained in the system allows evaluation of the 3qSS parameters. In particular the arrival rates of the aggregated class are estimated as an average of arrival rates of classes weighted with their steady state probabilities; in order to ensure convergence, the condition of equal throughput among the different synchronization stations is also used. See Appendix A for the details of the method. Thus, there is a total of R þ 1 synchronization stations that globally model the entire synchronization station. In particular there are R þ 1 3qSS modeling the behavior of class r customers (where r ¼ 1; . . . ; R) and one 2qSS modeling the behavior of class R þ 1. Note that for the last class it is possible to properly model its behavior using a 2qSS model [8] since arrival rates keD ðneD Þ depend only on queue ED ; therefore, it is possible to aggregate in one class the behavior of the other classes r (where r ¼ 1; . . . ; R). The 3qSS is solved numerically using the SSOR algorithm [20]. Starting from an initial solution, the arrival rates of each 3qSS are changed until their throughputs converge to the same unknown value. The solution algorithm is similar to that proposed by Baynat and Dallery in [3] and described in Section 4.2.1. Also in this case, the major benefit of the technique is its lack of complexity since it requires Ri solutions of

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

321

3qSS and one solution of a 2qSS instead of solving one synchronization station with Ri þ 1 queues. See Appendix A for more details on the proposed technique.

5. Numerical tests 5.1. Accuracy of the method Our goal in this subsection is to measure the accuracy of the applied methodology in evaluating the performance of assembly systems controlled with kanbans. To do this we test the proposed methodology in three different cases. In the first case, shown in Fig. 11, the behavior of a system made up of two phases is analyzed. The first stage contains an assembly station where two components are combined in a single product that, after having completed the second stage, is sold as a finished product. The number of kanbans is the same for both components: 3 in the first stage and 5 in the second stage. Stations are exponential single servers with unit service times and a FIFO queue management rule. We assume for simplicity that raw components are always available and therefore we eliminate queues FP0;1 and FP0;2 before stage S1 . Different values of customer demand rates are considered so that the demand is always less than the maximum throughput of the systems respectively equal to 0.849 ± 10
3 (independent kanban systems) and 0.818 ± 10
3 (simultaneous kanban systems) calculated with simulation. In order to evaluate the accuracy of the proposed technique, the numerical results of the analytical methods are compared with the simulation results. The performance indicators we consider are the average throughput of the system (number of parts per time unit), the average level of finished parts, the average number of customer orders delays and the average service level (percentage of order satisfaction). Both types of assembly systems, independent and simultaneous, are analyzed with the approximation technique for general closed queuing networks. To analyze independent synchronization stations in isolation we use the approximation method to calculate the performance of a multi-class synchronization station described in full details in Appendix A. Furthermore, to evaluate the accuracy of the approximate method for synchronization stations, we use the exact numerical solution of the complete Markov chain. We use the notation independent approximate (3qSS) and independent approximate (exact) respectively, for

k 1,1 D1,1

FP1,1 S 1,1

S2 FP2

Finished Products

S 1,2 D3

FP1,2

Demand

D1,2 k 1,2

D2 k2

Fig. 11. Case 1: a two-stage assembly system with two components.

322

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

referring to performance indicators of independent assembly systems calculated using the approximated technique for the synchronization station proposed in this paper or by solving it exactly. Correspondingly, notations simultaneous approximate (3qSS) are used to refer to performance indicators of simultaneous assembly systems calculated using the approximated technique proposed in this paper for analyzing synchronization stations in an isolation mode. Simultaneous approximate (exact) is used to refer to performance indicators of simultaneous assembly systems calculated by analyzing exactly synchronization stations in an isolation mode by application of traditional numerical techniques. Furthermore, the notation simultaneous approximate (2qSS) is adopted when the Baynat and Dallery’s technique [3] is used to calculate the performance of synchronization stations in an isolation mode. Fig. 12 shows, for the first case, the service level and the average time delay of customer orders for the independent system as a function of customer demand rate. In this particular case we apply only simulation and approximate analytical methods that use the exact solution for analyzing the assembly station. The reason is quite simple; the assembly station is made up of only three queues and, therefore, solving the assembly station with the decomposition technique would require solving of three different blocks made up

INDEPENDENT (SIMULATION) INDEPENDENT (EXACT SS)

Service Level

100% 90% 80% 70% 60% 50% 40% 30% 20% 10% 0% 0.2

0.3

0.4

0.5

0.6

0.7

0.75

0.8

0.81

0.82

0.83

0.84

Demand rate INDEPENDENT (SIMULATION) INDEPENDENT (EXACT SS)

Time Delay

90 80 70 60 50 40 30 20 10 0 0.2

0.3

0.4

0.5

0.6

0.7

0.75

0.8

0.81

0.82

0.83

0.84

Demand rate Fig. 12. Case 1: average service level and time delay of customer orders for the independent assembly system.

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

323

of three queues. Since each block is analyzed by solving exactly its Markov chain, it is convenient to directly solve the assembly station once. As the demand rate decreases, the service level approaches 100% since it is now easier for the system to satisfy customer orders. Results from analytical methods are very close to those for simulation, in particular in the left area of the graph in which the service level values are high. This area represents the most realistic cases in which customer orders must be satisfied with probabilities that are typically greater than 80%. In this particular case, the errors calculated on simulation results are small and less than 8% for the service level: however for small values of service levels, errors can be larger. Note that the accuracy of analytical methods increases when customer demand rate decreases. Fig. 13 shows the service level and the average time delay of customer orders for the simultaneous system as a function of customer demand rate. Also in this case it can be noticed that results are accurate in the realistic cases, that is when the system service level assumes values greater than 80%.

SIMULTANEOUS (SIMULATION) SIMULTANEOUS (EXACT SS)

Service Level

100% 90% 80% 70% 60% 50% 40% 30% 20% 10% 0% 0.2

0.3

0.4

0.5

0.6

0.7

0.75

0.8

0.81

0.82

0.83

0.84

Demand rate SIMULTANEOUS (SIMULATION) SIMULTANEOUS (EXACT SS)

Time Delay

100 90 80 70 60 50 40 30 20 10 0 0.2

0.3

0.4

0.5

0.6

0.7

0.75

0.8

0.81

Demand rate Fig. 13. Case 1: average service level and time delay of customer orders for the simultaneous assembly system.

324

Table 1 Case 1: results for demand rate equal to 0.2, 0.4 and 0.5 Model

Service level

FP1;1

FP1;2

FP2

D3

Time delay

0.2

Simultaneous simulation Simultaneous approximate (exact) Independent simulation Independent approximate (exact)

0:9997  10
4

2:7524  10
4

2:7525  10
4

4:7476  10
4

0

0:0004  10
4

0.9997

2.7525

2.7525

4.7463

8:6  10
5

0.0004

0:9997  10
4

2:7519  10
4

2:7521  10
4

4:7476  10
4

0

0:0004  10
4

0.9997

2.7520

2.7520

4.7463

8:6  10
5

0.0004

Simultaneous simulation Simultaneous approximate (exact) Independent simulation Independent approximate (exact)

0:9865  10
4

2:3850  1:4  10
4

2:3844  1:4  10
4

4:2855  10
4

0:0100  10
4

0:0250  1:1  10
4

0.9858

2.3834

2.3834

4.2629

0.0105

0.0262

0:9867  10
4

2:3776  10
4

2:3779  10
4

4:2904  10
4

0:0097  10
4

0:0242  10
4

0.9861

2.3763

2.3763

4.2663

0.0102

0.0256

Simultaneous simulation Simultaneous approximate (exact) Independent simulation Independent approximate (exact)

0:9516  7  10
5

2:1425  2  10
4

2:1425  2  10
4

3:8671  10
4

0:0585  2  10
4

0:1170  3  10
4

0.9485

2.1419

2.1419

3.8170

0.0625

0.1251

0:9540  5  10
5

2:1302  2  10
4

2:1296  2  10
4

3:8824  4  10
4

0:0533  10
4

0:1066  10
4

0.9505

2.1268

2.1268

3.8339

0.0589

0.1179

0.4

0.5

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

Demand rate

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

325

Detailed results on case 1 are reported in Table 1 for different values of customer demand. Computational time is in the order of minutes and analytical methods always perform faster than simulation. The second test case, shown in Fig. 14, considers an assembly system made up of two different stages. The first stage contains an assembly station where three components are combined into a single product that, after having been worked in the second stage, is sold as a finished product. Also in this case stations are exponential single servers with unit service times and the FIFO queue management rule. We assume for simplicity that raw components are always available and therefore we eliminate queues FP0;1 , FP0;2 and FP0;3 before stage 1. Demand rate is constant and equal to 0.75. However, different combinations of numbers of kanbans in the two stages are considered when system performance is evaluated. Numerical results shown in Table 2 confirm case 1. Analytical methods provide accurate results for large values of service levels (i.e. when there is a large number of kanbans in the system). Analytical methods perform generally faster than simulation. In cases where convergence is too long (e.g. large number of kanbans) the algorithm exits from the iterative procedure and results are taken from the last iteration; in this case also the analytical methods provide accurate results. A more complex system is considered in the third case (see Fig. 15) where we analyze an assembly system composed of three different stages. Note that the second stage of the system, crossed by components of type 1, is made up of different machines. Also in this case stations are exponential single servers with values of service rates reported in the figure and FIFO queue management rule. The number of kanbans is equal to 3 for each stage. Results from simulation and analytical methods are reported in Fig. 16 for different values of demand rate. The errors made by the analytical method increase in this last case due to the complexity of the queuing network. However, analytical methods provide fairly accurate results for large values of service levels corresponding to the most realistic cases. Notice also that the proposed analytical method seems to be more accurate than the one in [3], denoted with 2qSS, in analysis of simultaneous assembly systems. Table 3 reports the detailed results for two different values of demand rate. In the first scenario the demand rate is high in comparison with the maximum throughput of the systems equal to 0.1854 ± 10
4 and 0.1853 ± 10
4 for independent and simultaneous systems respectively calculated with simulation.

k 1,1

S1,1 D1,1

FP1,1 S2

S 1,2 D1,2

FP 2

FP1,2

Finished Products

k 1,2 FP1,3

D3

S1,3

Demand

D1,3 k 1,3 D2 k2

Fig. 14. Case 2: a two-stage assembly system with three components.

326

Table 2 Case 2: results for different numbers of kanbans Model

Service level

D1;1

FP1;1

FP1;2

FP1;3

FP2

D3

Time delay

10, 10, 10, 8

Simultaneous simulation Simultaneous approximate (3qSS) Simultaneous approximate (2qSS) Independent simulation Independent approximate (3qSS)

0:8724 1:5  10
4 0.8537

1:9693 7:1  10
4 2.7462

7:2807 7:5  10
4 7.2538

7:2807 7:5  10
4 7.2538

7:2807 7:5  10
4 7.2538

5:0587 7:5  10
4 4.8863

0:4471 1:2  10
3 0.5510

0:5961 1:5  10
3 0.7347

0.8284

2.8195

7.1804

7.1804

7.1804

4.7036

0.7298

0.9730

0:8750 10
4 0.8576

2:0152 4  10
4 2.7999

7:235 5  10
4 7.2001

7:2373 5  10
4 7.2001

7:2386 5  10
4 7.2001

5:0768 8  10
4 4.9161

0:4274 8  10
4 0.5255

0:5700 1:1  10
3 0.7007

Simultaneous simulation Simultaneous approximate (3qSS) Simultaneous approximate (2qSS) Independent simulation Independent approximate (3qSS)

0:9059 1  10
4 0.8765

1:8483 5  10
4 2.6085

5:4013 6  10
4 5.3915

5:4059 6  10
4 5.3915

5:4022 6  10
4 5.3915

6:6509 1  10
3 6.2904

0:3593 7  10
4 0.5481

0:4790 9  10
4 0.7308

0.8323

2.6936

5.3064

5.3064

5.3064

5.8707

0.9229

1.2305

0:9122 9  10
5 0.8910

1:9135 4  10
4 2.6826

5:3363 5  10
4 5.3174

5:3302 8  10
4 5.3174

5:3340 5  10
4 5.3174

6:7151 1  10
3 6.4403

0:3159 8  10
4 0.4465

0:4212 1  10
3 0.5953

Simultaneous simulation Simultaneous approximate (3qSS) Simultaneous approximate (2qSS) Independent simulation Independent approximate (3qSS)

0:9258 1  10
4 0.905845

1:9863 8  10
4 2.83678

7:2640 9  10
4 7.16322

7:2549 6  10
4 7.16322

7:2620 7  10
4 7.16322

6:9126 1  10
3 6.71707

0:2557 8  10
4 0.486628

0:3410 1  10
3 0.364971

0.8865

2.8285

7.1715

7.1715

7.1715

6.4621

0.4984

0.6645

0:9266 1  10
4 0.9124

2:0379 9  10
4 2.8131

7:2119 9  10
4 7.1870

7:2105 1  10
3 7.1870

7:2128 7  10
4 7.1870

6:92181  10
3

0:2510 9  10
4 0.3253

0:3350 7  10
4 0.4338

8, 8, 8, 10

10, 10, 10, 10

6.7294

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

ðk1;1; k1;2 ; k1;3 ; k2 Þ

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

327

k 2,1 k 1,1 µ=0.6

D1,1

D 2,1

µ=0.5

0.5

µ=0.2

0.5

µ=1

FP2,1

µ=0.6

FP1,1 µ=1 D1,2

FP1,2

µ=1

µ=0.9

Finished Products FP3

FP2,2 D2,2 k2,2

D4

k 3,2

Demand

µ=0.6

D2,3

FP2,3

k2,3 D3 k3

Fig. 15. Case 3: a three-stage assembly system with three components.

5.2. Independent vs. simultaneous Independent systems are more reactive to the behavior of market demand because they transmit more quickly, compared with simultaneous systems, information on demand consumption to upstream stages. Fig. 17 shows the percentage difference on the average time delay of customer orders and on the average service level between independent and simultaneous systems calculated both with simulation and analytical methods. In particular the values represented in the figure are calculated in the following way: Percentage difference ¼

Independent
Simultaneous  100: Simultaneous

The values of service levels for both systems are reported in Fig. 12 and Table 1. It can be noted that independent systems seem to perform better than simultaneous ones for high saturation of the system, i.e. when the system throughput assumes values close to the demand rate. Therefore, in this particular case, the choice of the right control rule for the release of kanbans to be used in assembly synchronization stations may be important.

6. Conclusions This paper presented assembly systems controlled with kanbans. In particular, two different kanban policies, differing in the way in which kanbans are released, have been considered: independent and simultaneous systems. Both systems, modeled by means of queuing networks, have been assessed and compared with approximate analytical methods.

328

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336 IN DEPEND ENT (SIM ULAT ION) IN DEPEND ENT APPROXIMAT E (3 qSS)

Service level

1.0000 0.9800 0.9600 0.9400 0.9200 0.9000 0.8800 0.05

0.06

0.07

0.08

0.09

0.1

Demand rate SIMU LT ANEOUS (SIMUL ATION ) SIMU LT ANEOUS APPR OXI MATE (3qSS) SIMU LT ANEOUS APPR OXI MATE (2qSS)

Service level

0.9900 0.9700 0.9500 0.9300 0.9100 0.8900 0.8700 0.8500 0.8300 0.8100

0.05

0.06

0.07

0.08

0.09

0.1

Demand rate Fig. 16. Case 3: average service level for independent and simultaneous assembly systems.

Results reported in the paper demonstrate the accurateness of the developed technique in some defined cases where the service levels of the system assume large values close to 100%. Furthermore, the two different control kanban rules are compared in the paper. Systems functioning with the two analyzed rules perform differently on a high-saturated demand level, which is when customer demand requests are close to the maximum throughput of the system. On the other hand, both rules provide similar results when service levels are very high, thus meaning that either rule can be selected in these cases. Appendix A. Analysis of independent synchronization stations in an isolation mode A.1. Decomposition of the system In this Section we describe in detail the approximate analytical method that we developed to calculate the performance of synchronization stations with independent release of kanbans. Section 4 has showed

Demand rate

Model

Service level

FP1;1

FP1;2

FP2;1

FP2;2

FP2;3

FP3

Time delay

0.1

Simultaneous Simulation Simultaneous approximate (3qSS) Simultaneous approximate (2qSS) Independent simulation Independent approximate

0:9674 1  10
4 0.9050

2:8909 1  10
4 2.8890

2:8914 1  10
4 2.8890

1:7527 1  10
4 1.5531

2:8781 1  10
4 2.8770

2:8070 1  10
4 2.8053

2:7071 1  10
4 2.4829

0:34711 1  10
4 1.7236

0.8241

2.8890

2.8890

1.5320

2.8751

2.8009

2.2291

5.1467

0:9688 1  10
4 0.9082

2:8910 1  10
4 2.8890

2:8895 1  10
4 2.8890

1:7492 1  10
4 1.5517

2:87591 10
4 2.8752

2:8012 1  10
4 2.8012

2:7071 1  10
4 2.4944

0:3351 1  10
4 1.6335

Simultaneous simulation Simultaneous approximate (3qSS) Simultaneous approximate (2qSS) Independent simulation Independent approximate (exact)

0:9992 1  10
4 0.9950

2:9477 1  10
4 2.9474

2:9477 1  10
4 2.9474

2:4796 1  10
4 2.3851

2:9417 1  10
4 2.9412

2:9092 1  10
4 2.9092

2:9370 1  10
4 2.9105

0:0060 1  10
4 0.0450

0.9948

2.9474

2.9474

2.3839

2.9412

2.9091

2.9097

0.0631

0:99931 1  10
4 0.9962

2:9475 1  10
4 2.9474

2:9476 1  10
4 2.9474

2:4757 1  10
4 2.3846

2:9414 1  10
4 2.9412

2:9098 1  10
4 2.9092

2:9376 1  10
4 2.9170

0:0037 1  10
4 0.0394

0.05

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

Table 3 Case 3: numerical results for values of demand rate equal to 0.1 and 0.06

329

330

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336 SIMULATION EXACT SS

Service Level

1.60% 1.40% 1.20% 1.00% 0.80% 0.60% 0.40% 0.20% 0.00%

0.2

0.3

SIMULATION

0.4 Demand rate

EXACT SS

0.5

0.6

Time Delay

0.00% -2.00%

0.2

0.3

0.4

0.5

0.6

-4.00% -6.00% -8.00% -10.00% -12.00% -14.00% -16.00% -18.00%

Demand rate Fig. 17. Case 1: independent vs. simultaneous.

how a synchronization station with Ri components can be transformed into an equivalent and simpler Markovian model. However, the equivalent model is still too complex to be solved numerically. Therefore, an extension of the class-aggregation technique, introduced by Baynat and Dallery for analyzing multi-class service stations [3], has been developed and described in the remainder of Appendix A. Let us consider the equivalent system represented in Fig. 18 and composed of R þ 1 queues (for simplicity’s sake, we omit, in this Section, the index i related to the original synchronization station Ji;iþ1 ) each one visited by a different class of customers. Arrival rates depend on the buffer levels of the system. In particular, parameters ker ðner ; neD Þ (where ner and neD are limited by the maximum number of kanbans, i.e. ner ¼ 0; 1; . . . ; Nre and neD ¼ 0; 1; . . . ; NDe ) depend on the levels of queues Er and ED while keD ðneD Þ depends only on the level of queue ED . Each class r (where r ¼ 1; . . . ; R; R þ 1) sees the other customers as a resource since at least one customer of each class must be present at the same time in order to proceed to the next station. The vector  ne ¼ ðne1 ; . . . ; neR ; neD Þ describes the possible states of the system. Since a synchronization operation is assumed to be instantaneous, at least one queue of the system is empty at any time t, i.e. at least one element of  ne is equal to 0.

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

λ1e (n1e , nDe )

e

λ …e (n…e , nD )

λ Re (nRe , nDe)

λ D (n De ) e

331

E1

E :::

ER

Assembled part and (i+1)-stage kanban ED

Fig. 18. Equivalent model of the independent assembly synchronization station Ji;iþ1 .

The goal of the analysis is to calculate the stationary probabilities P ¼ ðne1 ; . . . ; neR ; neD Þ, on the basis of which the throughput of the synchronization station can be evaluated for all classes. To do this we decompose (see Fig. 19) the synchronization station with R þ 1 queues into a set of simpler blocks that reproduce, in an approximate manner, the behavior of the original synchronization station. Blocks are composed of three queues: the first queue represents the behavior of class r in queue Er , the second queue represents that of class R þ 1 (i.e. the class related with the demand), while the third queue represents aggregately that of the other classes not considered by the first two queues (see Fig. 20). We indicate by nrr , nra and nrD the buffer levels of the 3qSS synchronization station related to class r, where krr ðnrr ; nrD Þ, kra ðnra ; nrD Þ and krD ðnrD Þ are their corresponding arrival rates. E11

Class 1 Ea1

Aggregated Class

Synchronized part

1 D

E

Class 4

R=3

E22

E1 Class 2

Class 1

Ea2

Aggregated Class

E2

Class 2

Synchronized part

2 D

E

Class 4

E3

Synchronized part

Class 3

3

E3

Class 3 3

Di+1

Aggregated Class

Class 4

Ea E

3 D

Synchronized part

Class 4

Demand

k i+1

ED

Class 4 D

Aggregated Class

Ea

Synchronized part

Fig. 19. Decomposition of a synchronization station with four classes.

332

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

r

r

λrr(nr , nD )

Er λar (nar , nDr )

λar(nar, nDr)

Ea λDr(nDr )

λDr (nDr)

Ea

ED

ED

For r=1,…, R

For r=R+1

Fig. 20. Examples of building blocks: 3qSS and 2qSS.

Similarly with the approach adopted by Baynat and Dallery [3], the parameters of blocks (i.e. arrival rates and queues capacities) must be selected appropriately so that throughput of parts leaving the station is equal for all customer classes. A set of equations concerning the equality of throughputs among the blocks allows evaluation of the block parameters. Thus, there is a total of R þ 1 synchronization stations that globally model the entire synchronization station. In particular there are R þ 1 3qSS modeling the behavior of class r customers (where r ¼ 1; . . . ; R) and one 2qSS modeling the behavior of (R þ 1)th class. Note that the last class can be properly modeled using a 2qSS model because arrival rates keD ðneD Þ depend only on queue Ed (see Fig. 20). Therefore, it is possible to aggregate in one class the behavior of the other classes r with r ¼ 1; . . . ; R. In this case, the block related to class R þ 1 is a synchronization station that has only two queues. The first queue models the behavior of class R þ 1 reproducing in the same way as queue ED of the equivalent system. The second queue models in an aggregate way the behavior of all the other classes s (where s ¼ 1; . . . ; R). We use the method adopted by Baynat and Dallery [3] for analyzing the 2qSS model.

A.2. Determination of block parameters In this Section we derive a set of equations that can be used to determine the unknown parameters of blocks that are derived from the decomposition of the system: arrival rates and maximum number of customers of queues: Er , Ea and ED : If the arrival rates into block r were known, it would be possible to solve numerically the Markov chain related to the single sub-system. Then the stationary probabilities P ðnrr ; nra ; nrD Þ could be evaluated, allowing calculation of the throughput Xrr of the class r as Xrr

¼

r N r
1 X

krr ðnrr ÞP r ðnrr Þ:

ðA:1Þ

nrr ¼0

In order to represent the original system properly, throughputs of each synchronization station must be equal to the same value; thus, the following relationships must be satisfied: Xrr ¼ X

with r ¼ 1; . . . ; R þ 1;

where X is the throughput of the station. The problem is limited to finding proper values for the arrival rates and the maximum number of customers in the blocks. The maximum number of customers in the block corresponding to class r ¼ 1; . . . ; R is Nrr ¼ Nre

with r ¼ 1; . . . ; R;

NDr

with r ¼ 1; . . . ; R:

¼

NDe

ðA:2Þ

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

333

The maximum number of parts that can be present in the queue containing customers of the aggregated class ar is limited by the lowest population size of aggregated classes and is given by Nar ¼

fNse g:

min

ðA:3Þ

s¼1;...;R s6¼r;s6¼Rþ1

The arrival process of class r in the rth block is the same as the original synchronization station:   krr ðnrr ; nrD Þ ¼ ker nrr ; nrD Þ ¼ ker ðnrr
minfnrr ; nrD g with nrr ¼ 0; 1; . . . ; Nre
1; nrD ¼ 0; 1; . . . ; NDe
1

and

r ¼ 1; . . . ; R:

ðA:4Þ

krr ðNre ; nrD Þ

Note that is equal to 0, that is no unit can arrive because the parts of class r are all contained in queue Er . In the same way the arrival process of class R þ 1 in the rth 3qSS is the same as the original synchronization station: krD ðnrD Þ ¼ keD ðnrD Þ

with nrD ¼ 0; 1; . . . ; NDe
1: kra ðnra ; nrD Þ

ðA:5Þ nra

1; . . . ; Nar

Now we have to evaluate only parameters for ¼
1; r ¼ 1; . . . ; R and r ¼ R þ 1. Let  nr ¼ ½ns  (with s ¼ 1; . . . ; R and s 6¼ r, s 6¼ R þ 1) denote the state-vector that represents the number of customers of all the classes except types r and R þ 1, and with Xðnra ; nrD Þ the set of all possible vectors nr with the same associated number of class ar and R þ 1 vectors: n o Xr ðnra ; nrD Þ ¼  nr ¼ ½ns s¼1;...;R such that ns P nra and 9t; t 6¼ r and t 6¼ R þ 1 such that nt ¼ nra : s6¼r

The set contains all the possible combinations ðnra ; nrD Þ that make possible a number nra of parts in the aggregated queue of the 3qSS and a number nrD in the queue of class R þ 1. Using this definition, the arrival rates of class ar customers corresponding to a particular state ðnra ; nrD Þ are as follows:  e ks ðns ; nrD Þ if ns ¼ nra and nt > nra for all t 6¼ s and t 6¼ r and t 6¼ D; kra ð ðA:6Þ nr ; nrD Þ ¼ 0 otherwise: If all the original queues nt are greater than nra except for queue ns (i.e. only one unit of class s is required to allow synchronization of parts in the aggregated queues), then the arrival rate is equal to ks ðns Þ. If at least two queues are equal to nra , then synchronization of parts in the aggregated queues is not possible. Thus, the values kra ðnra ; nrD Þ can be calculated as a weighted average of arrival rates kra ðnr ; nrD Þ: P kra ð nr ; nrD ÞP r ð nr ; nrD Þ kra ðnra ; nrD Þ ¼

ð nr ;nrD Þ2Xr ðnra ;nrD Þ

P

ð nr ;nrD Þ2Xr ðnra ;nrD Þ

P r ð nr ; nrD Þ

ðA:7Þ

;

where the weights correspond to the probabilities of being in that state. The sum of all the probabilities P r ð nr Þ such that nr 2 Xr ðnra ; nrD Þ can be expressed as 0 1 B C B r C P r ð nr ; nrD Þ ¼ P r B n ¼ ½ns  s6¼r such that ns P nra ; nrD C @ A s6¼Rþ1 ð nr ;nr Þ2Xr ðnr ;nr Þ X

D

a

D

0

1

B C B r C
P r B n ¼ ½ns  s6¼r such that ns > nra ; nrD C: @ A s6¼Rþ1

ðA:8Þ

334

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

Substituting (A.6) and (A.8) in Eq. (A.7), the arrival rates kra ðnra ; nrD Þ can be calculated as 0 P ð nr ;nrD Þ2Xr ðnra ;nrD Þ

kra ðnra ; nrD Þ ¼

B r kra ð nr ; nrD ÞP r @ n ¼ ½ns 

s6¼r s6¼Rþ1

0 B r P r @ n ¼ ½ns 

C : ns ¼ nra and nt > nar and t 6¼ r and s; nrD A

1 s6¼r s6¼Rþ1

1

0

1

C B : ns P nra ; nrD A
P r @nr ¼ ½ns 

s6¼r s6¼Rþ1

:

C : ns > nra ; nrD A

To calculate the stationary probabilities P r we introduce the approximation   Y Pr  Pss ðns ÞP D ðnrD Þ; nr ¼ ½ns ; nrD ¼ s6¼r s6¼Rþ1

where Pss ðns Þ represents the marginal probability that ns customers of class s are present in the synchronization station. Since these probabilities are the unknown quantities of the problem, i.e. the probabilities that we want to evaluate, it is necessary to adopt an iterative approach. In particular the values of the last iteration are used as an estimation of Pss ðns Þ. In such a way we can write Q P s r r ks ðna ; nD Þ  Pss ðns ¼ nra Þ  P D ðnrD Þ  Ptt ðnt nra Þ kra ðnra ; nrD Þ

¼

s6¼r s6¼Rþ1

Q

s6¼r s6¼Rþ1

Pss ðns nra Þ  P D ðnrD Þ


Q s6¼r s6¼Rþ1

t6¼s t6¼r;Rþ1

ðA:9Þ

Pss ðns > nra Þ  P D ðnrD Þ

and obviously kra ðNar ; nrD Þ ¼ 0. Likewise, it is possible to extend Eqs. (A.2) and (A.7) to calculate the arrival rates of the last block corresponding to class R þ 1. Relation (A.9) does not work well when nra ¼ 0 due to the approximation in the estimation of the detailed probabilities P r ðnra ; nrD Þ by the product of the marginal probabilities. This approximation is less accurate for low values of nra . Thus, to determine the parameters kra ð0; nrD Þ we oblige the throughputs Xrr (to be calculated from relation (A.1)) to be equal by explicitly using the following equations to calculate the arrival rates of the aggregate classes: Xrr ¼ X

with r ¼ 1; . . . ; R þ 1:

Furthermore, the throughput can be calculated with the following equation: X ¼

Rþ1 X   krr ð0ÞP r nrr ¼ 0; nra > 0; nrD > 0 : r¼0

Finally we have

PRþ1

r r¼1 Na

for r ¼ 1; . . . ; R þ 1

ðA:10Þ

equations to determine the same number of arrival rates of aggregated classes: P s r r s Q t r D r r

8 > > > < > > > :

ks ðna ;nD ÞPs ðns ¼na ÞP ðnD Þ

kra ðnra ; nrD Þ ¼

s6¼r s6¼Rþ1

Q

s6¼r s6¼Rþ1

Pss ðns nra ÞP D ðnrD Þ


Xrr ¼ X

Q

Pt ðnt P na Þ

t6¼s t6¼r;Rþ1 Pss ðns >nra ÞP D ðnrD Þ

:

ðA:11Þ

s6¼r s6¼Rþ1

The method has been implemented as follows. At step 0 all the parameters of the algorithm are initialized; then each decomposed block is analyzed at step 1 by applying the SSOR algorithm [20]. At step 2 the throughput of the station is calculated and the exit condition is checked at step 3; the algorithm loops again to step 1 if the exit condition is not satisfied.

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

335

The algorithm Step 0 Initialization Calculate Nrr Nar and NDr from relations (A.2) and (A.3). Calculate krr ðnrr ; nrD Þ and krD ðnrD Þ from relations (A.4) and (A.5). Assign initial values to the following coefficients: k ¼ 0, where k represents the number of iterations CðkÞ ¼ 1, where CðkÞ represents a convergence coefficient at the iteration k. Step 1 Analysis of blocks k ¼kþ1 For r ¼ 1; . . . ; R þ 1 (a) If k ¼ 1 then initialize kra ð0; nrD Þ from relation (A.9) (b) Else calculate kra ð0; nrD Þ from the following relation: 

 1 Xavg ðkÞ 1 Xavg ðk
1Þ r k ð0; nrD ÞðkÞ kra ðnra ; nrD Þ ¼ þ 1
 r CðkÞ X ðk
1Þ CðkÞ Xr ðk
1Þ a where Xavg ðkÞ is the average throughput at iteration k: Xavg ðkÞ ¼

Rþ1 1 X X r ðkÞ R þ 1 r¼1 r

(c) Solve the Markov chain of the rth block (d) Calculate the marginal probabilities Prr ðnrr ; nrD Þ (e) Calculate the class-r throughput Xrr ðkÞ Step 2 Calculation of throughput (a) Calculate the throughput X ðkÞ from relation (A.10) (b) Calculate the average throughput Step 3 Exit condition (a) Calculate the coefficient precision as:  X ðkÞ
X ðk
1Þ X ðkÞ
Xrr ðkÞ ; precisionðkÞ ¼ max r¼1;...;Rþ1 X ðkÞ X ðkÞ (b) If precisionðkÞ < e then exit. Else If precisionðkÞ > precision ðk
1Þ then CðkÞ ¼ CðkÞ þ 1 If precisionðkÞ > e and k > 1 then 

precisionðkÞ
precisionðk
1Þ CðkÞ ¼ max Cðk
1Þ þ 100  precisionðk
1Þ Goto step 1.

References [1] B. Baynat, Y. Dallery, Approximate techniques for general closed queueing networks with subnetworks having population constraints, European Journal of Operational Research 69 (1993) 250–264.

336

A. Matta et al. / European Journal of Operational Research 166 (2005) 310–336

[2] B. Baynat, Y. Dallery, A unified view of product-form approximation techniques for general closed queueing networks, Performance Evaluation 18 (3) (1993) 205–224. [3] B. Baynat, Y. Dallery, Approximate analysis of multi-class synchronized closed queueing networks, in: MASCOTS’95, Proceedings of the Third International Workshop on Modeling, Analysis, and Simulation of Computer and Telecommunication Systems, 18–20 January 1995, pp. 23–27. [4] B. Baynat, Y. Dallery, M. Di Mascolo, Y. Frein, A multiclass approximation technique for the analysis of kanban-like control systems, International Journal of Production Research 39 (2) (2001) 307–328. [5] J.-M. Bollon, M. Di Mascolo, Y. Frein, Unified framework for describing and comparing the dynamics of pull control policies, Annals of Operations Research 125 (2004) 21–45. [6] J.A. Buzacott, I. Cheinis, PAC controlled assembly systems, in: Conference on Multi-Echelon Inventory Systems, IBM T.J. Watson Research Center, June 7–9, 1993. [7] J.A. Buzacott, J.G. Shanthikumar, Stochastic Models of Manufacturing Systems, Prentice-Hall, 1993. [8] Y. Dallery, X.R. Cao, Operational analysis of closed stochastic queueing networks, Performance Evaluation 14 (1) (1992) 43–61. [9] S.L. De Araujo, Sur l’analyse et le dimensionnement de syst
emes de production geres en kanban, Ph.D. Thesis, Laboratoire d’Automatique de Grenoble (LAG), Institut National Polytechnique de Grenoble (INPG), September 1993. [10] M. Di Mascolo, Y. Dallery, Performance evaluation of kanban controlled assembly systems, in: Symposium on Discrete Events and Manufacturing Systems of the Multiconference IMACS-IEEE/SMC CESA’96, Lille, France, July 1996. [11] M. Di Mascolo, Y. Frein, Y. Dallery, An analytical method for performance evaluation of kanban controlled production systems, Operations Research 44 (1) (1996) 50–64. [12] M. Di Mascolo, Y. Frein, B. Baynat, Y. Dallery, Queueing network modeling and analysis of generalized kanban systems, European Control Conference (ECC’93), Groningen, The Netherlands, 28 June–1 July, 1993. [13] C. Duri, Y. Frein, M. Di Mascolo, Comparison among three pull control policies: Kanban, base stock, and generalized kanban, Annals of Operations Research 93 (2000) 41–69. [14] Y. Frein, M. Di Mascolo, Y. Dallery, On the design of generalized kanban control systems, International Journal of Operations & Production Management 15 (9) (1995) 158–184. [15] W.C. Gordon, G.F. Newell, Closed queueing networks with exponential servers, Operations Research 15 (1967) 252–267. [16] F. Karaesmen, Y. Dallery, Performance comparison of pull type control mechanisms for multi-stage manufacturing, International Journal of Production Economics 68 (2000) 59–71. [17] G. Liberopoulos, Y. Dallery, A unified framework for pull control mechanism in multi-stage manufacturing systems, Annals of Operations Research 93 (2000) 325–355. [18] D. Mitra, I. Mitrani, Analysis of a novel discipline for cell coordination in production lines: Part I, Management Science 36 (1990) 1548–1566. [19] D. Mitra, I. Mitrani, Analysis of a novel discipline for cell coordination in production lines: Part II––stochastic demands, Operations Research 39 (1991) 807–823. [20] B. Stewart, Introduction to the Numerical Solution of Markov Chains, Princeton University Press, 1994.