# Modeling and Analysis of Reentrant Manufacturing Systems: Micro- and Macroperspectives

Authors: Fenglan He; Ming Dong; Dong Yang
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

In order to obtain the better analysis of the multiple reentrant manufacturing systems (MRMSs), their modeling and analysis from both micro- and macroperspectives are considered. First, this paper presents the discrete event simulation models for MRMS and the corresponding algorithms are developed. In order to describe MRMS more accurately, then a modified continuum model is proposed. This continuum model takes into account the re-entrant degree of products, and its effectiveness is verified through numerical experiments. Finally, based on the discrete event simulation and the modified continuum models, a numerical example is used to analyze the MRMS. The changes in the WIP levels and outflux are also analyzed in details for multiple re-entrant supply chain networks. Meanwhile, some interesting observations are discussed.


---

## Body

## 1. Introduction

In recent years,  factories and production systems have become larger and more complicated. The reentrant manufacturing system is a typical example, in which the work in process (WIP) repeatedly passes through the same workstation at different stages of the process routes. Figure1 gives the structure of a reentrant manufacturing process. Here, B represents the storage areas of each work center. In the large factories, no experiments can be done involving whole supply chains. Therefore, simulation models are developed, which can be used to substitute for the real systems. Especially the multiscale simulations of production flows in manufacturing systems have become a very important research topic recently. Scaling also plays a role in biomedical engineering field [1]. The production flows naturally show two time scales: a shorter one that describes single items processed through the individual workstations and a longer one that describes the response time of the whole factory [2]. In fact, two scales of flows, small and large, are well explained in fractals [3]. Furthermore, the representations of their bounds are also provided [4].Figure 1
The structure of a reentrant manufacturing system.Currently, most of the semiconductor manufacturing systems are described by the discrete models, the advantages of the discrete models are (1) some of the actual static problems can be directly reflected in the discrete models. (2) In reality, the methods of data collection are discrete. (3) There are many existing methods to solve the discrete models. Based on the theory of the discrete event dynamic systems, there are several discrete methods used to model the multiple reentrant semiconductor production flows.First, queuing networks can intuitively describe the discrete production process of the wafer production line and obtain the analytical expression of the performance evaluation. A methodology for supply chain inventory analysis and optimization was presented by linking production authorization (PA) strategy to queuing models [5]. The statistics of arrival flows of the fluid models in queuing systems has been studied [6, 7]. Dong and Chen [8] proposed a network of inventory-queue models for the performance modeling and analysis of an integrated supply chain network. For a single multiple semiconductor manufacturing system, S. Kumar and P. R. Kumar [9] focused on the analysis of queuing theory of the reentrant manufacturing systems. However, with the emergence of new modes of production, especially for the development of the large quantity and multi-class order production mode, modeling and analysis of the traditional queuing network models becomes more difficult, a simplified model of wafer manufacturing systems becomes less practical. On the other hand,  it is very difficult to solve the large-scale queuing network models. Most of the queuing models can only be used to evaluate the stability of some scheduling policies, while they are difficult to be used directly for the real-world supply chains. The queuing models are only used to evaluate the stability of some scheduling policies and can not be directly used for the actual production systems [10].Second, Petri net (PN) model has been widely used in modeling, simulating, analyzing, and controlling the discrete event dynamic systems. Compared with some other description tools, PN model is especially easy to describe concurrent phenomena and simulate the parallel systems. However, with the increasing complexity and size of the manufacturing systems, the complexity of the PN model analysis is also a corresponding increase. Lin et al. [11] established a model of the reentrant semiconductor production lines using Petri-nets and studied on the stability of the system using buffer-boundless approach. Dong and Chen [12] developed a modular modeling approach based on object-oriented predicate/transition nets (OPTNs) for the analysis of supply chain process models. Liu et al. [13] proposed a model of semiconductor manufacturing systems on the basis of object-oriented colored time Petri nets and proved the validity of the model by simulating different scheduling rules.Third, fluid networks model comes from traffic theory and was introduced by Newell [14] to approximately solve queuing problems. This method treats the discrete manufacturing process as the dynamic allocation of processing capacity of equipment. The model can greatly reduce the dimensions of the system, thereby can reduce the difficulty of analysis, while it has good robustness and can be applied to random environment. However, the expression of the model is obtained by the constraint equations, so it is difficult to apply to multi-class and high volume make-to-order production mode. Dai and Weiss [15] analyzed the relationship between the stability of the fluid model and the stability of the scheduling policies for the related queuing networks. Here, the stability of the fluid models is expressed by the boundedness of the fluid variables for a fixed influx. Göttlich [16] deduced the conservation laws under the form of ordinary differential and proved the existence of its solution based on the fluid models of supply chain networks.On the other hand, the continuous models also have several advantages in description of the production lines. That is, they are scalable, more detailed results that can be found as compared to fluid models, and more important, they are amenable to optimization and control. Recently, the continuous models have been applied to many fields and have achieved some significant research results. Anderson [17] established the basic continuous model for supply chains and described production flow of the systems using rate equations macroscopically. Lee et al. [18] studied on the supply chain simulation with discrete-continuous combined modeling approach. This could combine the wide applicability of the discrete event simulation (DES) and fast computation of the continuum models together. Armburster and Ringhofer [19] introduced the concept of materials’ density and established the continuous model of large-scale reentrant manufacturing systems and proposed new state equations so that the continuous model could be used to the real production systems. Compared with a discrete model, van den Berg et al. [20] verified the validity of the continuous model of simple serial production systems and then solved the optimal control problems with consideration of the demand growth.The continuous models can more accurately reflect the actual situation, while the corresponding discrete models need to be discretized in time or space, then takes a constant value in each discrete node to represent the state of a period of time. This method is the approximation of the real situations and can not fully reflect the actual situations. Using the common ground and similarities of the semiconductor manufacturing systems, the large-scale complex network can be decomposed into relatively small-scale simple problems. The multiscale methods can be more accurate modeling and analysis of the manufacturing systems. Armbruster and Ringhofer [21] extended the standard stochastic models by introducing the concept of a random phase velocity. This leads to the concepts of temperature and diffusion in the corresponding kinetic and fluid models for supply chains. Zou et al. [22] demonstrated certain features of equation free coarse-grained computation for a reentrant supply-chain model. They took the advantage of the time scale separation to directly solve coarse-grained equations through an equation-free computational approach. Unver et al. [23, 24] presented a continuum (traffic flow like) model for the flow of products through complex production networks, based on statistical information obtained from extensive observations of the system. The resulting model consists of a system of hyperbolic conservation laws, which exhibit the correct diffusive properties given by the variance of the observed data.In this paper, in order for better analysis of the multiple reentrant manufacturing systems, their modeling and analysis from both micro- and macroperspectives are considered. First, the discrete event simulation models and their basic algorithm are proposed. In order to describe the multiple reentrant semiconductor manufacturing systems more precisely, a modified continuous model is proposed, which can reflect how the reentrant degree of a product impacts on the system performance. Finally, the changes of the WIP levels and outflux are analyzed based on the discrete event simulation models and the modified continuous models.The structure of this paper is organized as follows: Section2 presents the discrete event simulation models and their basic algorithm. In Section 3, in order to describe the manufacturing systems more accurately, a modified continuous model is proposed, which introduces the concept of reentrant factor, and the effectiveness is also verified through a numerical example. In Section 4, the changes of the WIP levels and outflux are analyzed between the discrete event simulation models and the modified continuous models. Meanwhile, some interesting observations are discussed. Finally,  some conclusions and suggested further areas of investigation are given in Section 5.


## 2. Simulation Models

### 2.1. Discrete Event Simulation Models

In the discrete event simulation (DES) models, each item to be produced is treated as individual. The production process consists ofM stages (s1,s2,…,sM), anm denotes the time at which the lot number n arrives at stage sm, and enm denotes the time at which the lot number n exits at stage sm (or the lot number n arrivers at stage sm+1). τnm denotes the processing time that the lot number n takes at stage sm, which is known as the throughput time (TPT). The relationship between the arrival and exit times is given via the law(2.1)(a)enm=anm+τnm,(b)dP{τnm=r}=Tm(r,anm)dr,(c)anm+1=enm.
Here, Tm(r,a) denotes the time-dependent distribution of throughput times of stage sm, and P denotes the probability distribution of the processing time. The TPT is determined by the state of the supplier sm at time anm when lot number n arrives. The throughput time distribution Tm is usually dependent on the total number lots Wm(t) processed at stage sm at time t, the so-called work in progress (WIP). The WIP Wm(t) is computed as (2.2)Wm(t)=∑nH(t-anm)-H(t-anm+1),
where H denotes the usual Heaviside function.In engineering practices, the most significant influence on the form of this distribution comes from the total number of items in progress, that is, the work in progress (WIP). So the TPT distributionTm in (2.1) is written as Tm=Tm(r,Wm(anm)). When the TPT at the beginning of the simulation is fixed, we essentially treat the whole factory as a single queue whose length at the item arrival time determines the time that the item needs to get processed. In the case of complicated supply chain networks, where each individual supplier models a whole factory, the dependence of the distribution Tm on the WIP Wm is more complex and could be given by experimental data.For a high-volume multiple reentrant manufacturing system, the time interval between jobs becomes less important, and then the jobs can be seen as a continuous way. In these continuous production models, the WIPWm(t) and the influx λm(t) are treated as a continuous function. This leads to a class of models often referred to as fluid models in the supply chain literature [25]. The primary variables are the WIP Wm(t) and the fluxes λm(t) from stage sm-1 to sm. A random throughput time function τ(t) can be computed from the WIP obtained before. The WIP Wm(t) is then given by the conservation law(2.3)ddtWm(t)=λm(t)-μm(t),λm+1(t)=μm(t).
Here, λ and μ represent the influx and outflux of the lots, respectively. According to (2.1), the fluxes are given (2.4)μm(t)=∫δ(s+τ(s)-t)λm(s)ds.For a large number of lots, the individual lots are replaced by a continuous quantity, and the fluxesλm are replaced by continuous functions. The WIP Wm(t) is evolved by the conservation law (d/dt)W(t)=λ(t)-μ(t). Using the discrete time steps, the following formulation can be obtained: (2.5)W(t+Δt)=W(t)+Δt[λ(t)-μ(t)].So the question can be used to solve outfluxμ(t) in terms of influx λ(t).However, if suppliersm is a reentrant system itself, then the lots arriving after time anm will significantly influence the TPT τnm since they will compete with lot number n in the individual subqueues of the system. The set of policies, such as FIFO, PUSH, or PULL, will be used to govern the competition. This situation is treated by introducing the concepts of phase (a scaled position) and phase velocity.In phase models, the evolution of a large ensemble of lots is modeled by describing the trajectories of each individual lot in phase space based on Newton equations. The trajectoriesξ(t) and the throughput time τ(t) are introduced. The throughput time τ(t) is a random variable sampled from the distribution T(r,t). Given a prescribed TPT distribution function T(r,t), which is expressed in terms of WIP as τ(r,WIP)), the position and the throughput time after a small enough time interval Δt are given by(2.6a)ξ(t+Δt)=ξ(t)+Δtτ(t),τ(t+Δt)=(1-k)τ(t)+kη(t),(2.6b)P{k(t)=1}=wΔt,P{k(t)=0}=1-wΔt,dP{η(t)=r}=T(r,t)dr,(2.6c)ξ(0)=0,dP{η(0)=r}=T(r,0)dr.A random variablek is used to decide whether to update the throughput time or not, which is either one or zero. The probability of k=1 equals wΔt and the probability of k=0 is 1-wΔt. Here, w is the update frequency. If k=1, the throughput time is updated and the new throughput time is given by a random number η(t), which is obtained from the throughput time distribution τ(r,WIP). If k=0, the current throughput time remains unchanged. In general, suppose that w depends on both η(t) and time t, then ((2.6b)) will be replaced by(2.6b′)P{k(t)=1}=Δtw(η(t),t),P{k(t)=0}=1-Δtw(η(t),t),dP{η(t)=r}=T(r,t)dr.
According to [18], w can be chosen as (2.7)w(r,t)=λ(t)rT-1,
where λ(t) is the influx of items and T-1=∫r-1T(r,t)dr. Based on the evolution (2.6a), (2.6b), and (2.6c), the algorithm procedure to solve the phase and TPT of an item for the discrete event and continuous production is then as follows.Step 1.
Set the initial phase to 0 as the item enters the factory, increase the WIPW(t) by one and adjust the probability distribution T(r,WIP) accordingly.Step 2.
Computew(r,t) based on (2.7), and sample the parameters k(t) and η(t) according to their distribution ((2.6b)).Step 3.
Compute the TPTτ(t+Δt) and the phase ξ(t+Δt) according to ((2.6b)).Step 4.
Update the distribution TPTT(r,WIP), according to the current WIP at the time t+Δt, and back to Step 2, repeat this loop until the phase ξ=1. Now the time when ξ=1 becomes the time of the item leaving the system.The functionξ(t) plays a role of a position which is advanced with a randomly changing velocity 1/η(t). It is, therefore, to derive a kinetic formulation. In the context of supply chains, these models are often referred to as traffic flow models [26]. A product arrives at the first stage at time t=a and moves from x=0 to x=1 along a trajectory with velocity 1/ηa(t), in such a way that it reaches the end x=1 at time t=e, then the velocity has to be chosen in the following way ∫ae(1/ηa(t))dt=1. The total WIP W(t) and the flux F(x,t) at any point x∈[0,1] are then obtained as follows: (2.8)W(t)=∫[H(ξa(t))-H(ξa(t)-1)]λ(a)da,F(x,t)=∫ηa(t)δ(x-ξa(t))λ(a)da,
where λ(a) is concentrated on the arrival times an for a discrete model or is a continuous function for a continuous production model.


### 2.2. Basic Partial Differential Equation (PDE) Models

Partial differential equation (PDE) models are actually continuum approximation of fluid models. Recently, PDE models for large-scale multiple reentrant production systems have become an important research topic. PDE models do have several advantages, that is, they might be preferred when evaluating the overall performance of large-scale manufacturing systems or solving the optimization problem for a longer time plan. Currently, the continuous models have been applied to many fields and have achieved some significant research results.This is appropriate to describe a semiconductor manufacturing fab involving a large number of items in many stages:ρ(x,t) is the density of the products with units (units/space) in the system, which is the conserved variable. Here, x denotes the degree of completion (DOC), x=0 describes raw products that have just entered into the factory, and x=1 denotes finished products that are ready to exit from the system. So, x is in the closed interval: x∈[0,1]. The total number of products in the system can be obtained by taking the integral of density ρ(x,t) of products over the stage variable x from 0 to 1. Therefore, the total WIP W(t) as a function of time can be obtained as follows: (2.9)W(t)=∫01ρ(x,t)dx.Assuming that there is a unique entry and exit for the system and the yield is 100%, PDE models can be given bellow according to the conservation law,(2.10)∂ρ(x,t)∂t+∂(v(ρ(x,t))ρ(x,t))∂x=0,x∈[0,1],t∈(0,∞),
where v(ρ(x,t)) is a velocity function that depends on the density ρ(x,t) only. For a reentrant supply chain system, we assume that v(ρ(x,t)) can be described by a state equation of the following form: (2.11)v(ρ(x,t))=v0(1-W(t)Wmax⁡),
where v0 is the velocity for the empty supply chain system, and Wmax is the maximal load (capacity of the supply chain system). Clearly, the velocity v(ρ(x,t)) is determined by the total WIP. The boundary condition for the start rate λ(t) of products entering the supply chain system at x=0 is then defined as(2.12)λ(t)=ρ(0,t)v(t).An arbitrary initial condition for the density of the products can be expressed as(2.13)ρ(x,0)=ρ0(x).The production process of the supply chain network can be described as an equivalentM/M/1 queue. The state variable ρeq denotes the equilibrium density of the supply chain system as a whole. Let ρ=ρeq, then ρ(0,t)=ρeq and v=veq=1/τ. Correspondingly, the associated cycle time in steady state is τ=1/veq. Since a job arriving at a queue with a processing rate is μ=1/vmax, the cycle time τ=μ(1+L) can be obtained according to the queuing theory, the equilibrium velocity therefore becomes(2.14)v(ρ)=vmax⁡1+∫01ρ(x,t)dx=vmax⁡1+w(t).Equation (2.14) is a widely used expression between v and ρ for large-scale multiple reentrant supply chain systems, notice that the velocity v only depends upon the WIP at stage x. Therefore, (2.10) can be rewritten as(2.15)∂ρ(x,t)∂t+v∂ρ(x,t)∂x=0,x∈[0,1],t∈(0,∞).Hence, assume that an initial WIP distributionρ0(x) in the factory is given, the resulting full PDE models for the single-product multiple reentrant supply chain systems are given as follows:(2.16)∂ρ(x,t)∂t+v(ρ)∂ρ(x,t)∂x=0,ρ(x,0)=ρ0(x),ρ(0,t)v(t)=λ(t),v(ρ)=vmax⁡1+∫01ρ(x,t)dx.If the influxλ(t) and the initial condition ρ0(x) are nonnegative, then the density will remain nonnegative. We use an upwinding scheme [27] to discretize the PDE which is given by the following equation:(2.17)ρ(xi,tj+1)=ρ(xi,tj)-ΔtΔxv(tj)[ρ(xi,tj)-ρ(xi-1,tj)],
where i=1,2,…,N, and j=0,1,…,M-1. Δt and Δx are the step sizes in time and space, respectively. Based on the boundary condition, the propagation scheme is given by(2.18)ρ(x0,tj+1)=ρ(x0,tj)-ΔtΔx[v(tj)ρ(x0,tj)-λ(tj)].Since the Courant-Friedrich-Levy (CFL) condition is necessary for stability, the time stepΔt and the space step Δx must satisfy the following formulation: |(Δt/Δx)vmax(t)|≤1. Here, vmax(t) is the maximum velocities in the system at time t.Based on the above formula, the density distribution of each momentρ(xi,tj) can be obtained. Therefore, the system throughput rate q(xN,tj) of each moment can also be computed as follows: (2.19)q(xN,tj)=ρ(xN,tj)v(ρ(xN,tj)).Furthermore, the total WIPw(t)=∫01ρ(x,t)dx can be obtained via the extended Simpson’s rule quadrature. Finally, the density distribution ρ(xi,tj) and throughput q(xN,tj) of each moment can be obtained.


## 2.1. Discrete Event Simulation Models

In the discrete event simulation (DES) models, each item to be produced is treated as individual. The production process consists ofM stages (s1,s2,…,sM), anm denotes the time at which the lot number n arrives at stage sm, and enm denotes the time at which the lot number n exits at stage sm (or the lot number n arrivers at stage sm+1). τnm denotes the processing time that the lot number n takes at stage sm, which is known as the throughput time (TPT). The relationship between the arrival and exit times is given via the law(2.1)(a)enm=anm+τnm,(b)dP{τnm=r}=Tm(r,anm)dr,(c)anm+1=enm.
Here, Tm(r,a) denotes the time-dependent distribution of throughput times of stage sm, and P denotes the probability distribution of the processing time. The TPT is determined by the state of the supplier sm at time anm when lot number n arrives. The throughput time distribution Tm is usually dependent on the total number lots Wm(t) processed at stage sm at time t, the so-called work in progress (WIP). The WIP Wm(t) is computed as (2.2)Wm(t)=∑nH(t-anm)-H(t-anm+1),
where H denotes the usual Heaviside function.In engineering practices, the most significant influence on the form of this distribution comes from the total number of items in progress, that is, the work in progress (WIP). So the TPT distributionTm in (2.1) is written as Tm=Tm(r,Wm(anm)). When the TPT at the beginning of the simulation is fixed, we essentially treat the whole factory as a single queue whose length at the item arrival time determines the time that the item needs to get processed. In the case of complicated supply chain networks, where each individual supplier models a whole factory, the dependence of the distribution Tm on the WIP Wm is more complex and could be given by experimental data.For a high-volume multiple reentrant manufacturing system, the time interval between jobs becomes less important, and then the jobs can be seen as a continuous way. In these continuous production models, the WIPWm(t) and the influx λm(t) are treated as a continuous function. This leads to a class of models often referred to as fluid models in the supply chain literature [25]. The primary variables are the WIP Wm(t) and the fluxes λm(t) from stage sm-1 to sm. A random throughput time function τ(t) can be computed from the WIP obtained before. The WIP Wm(t) is then given by the conservation law(2.3)ddtWm(t)=λm(t)-μm(t),λm+1(t)=μm(t).
Here, λ and μ represent the influx and outflux of the lots, respectively. According to (2.1), the fluxes are given (2.4)μm(t)=∫δ(s+τ(s)-t)λm(s)ds.For a large number of lots, the individual lots are replaced by a continuous quantity, and the fluxesλm are replaced by continuous functions. The WIP Wm(t) is evolved by the conservation law (d/dt)W(t)=λ(t)-μ(t). Using the discrete time steps, the following formulation can be obtained: (2.5)W(t+Δt)=W(t)+Δt[λ(t)-μ(t)].So the question can be used to solve outfluxμ(t) in terms of influx λ(t).However, if suppliersm is a reentrant system itself, then the lots arriving after time anm will significantly influence the TPT τnm since they will compete with lot number n in the individual subqueues of the system. The set of policies, such as FIFO, PUSH, or PULL, will be used to govern the competition. This situation is treated by introducing the concepts of phase (a scaled position) and phase velocity.In phase models, the evolution of a large ensemble of lots is modeled by describing the trajectories of each individual lot in phase space based on Newton equations. The trajectoriesξ(t) and the throughput time τ(t) are introduced. The throughput time τ(t) is a random variable sampled from the distribution T(r,t). Given a prescribed TPT distribution function T(r,t), which is expressed in terms of WIP as τ(r,WIP)), the position and the throughput time after a small enough time interval Δt are given by(2.6a)ξ(t+Δt)=ξ(t)+Δtτ(t),τ(t+Δt)=(1-k)τ(t)+kη(t),(2.6b)P{k(t)=1}=wΔt,P{k(t)=0}=1-wΔt,dP{η(t)=r}=T(r,t)dr,(2.6c)ξ(0)=0,dP{η(0)=r}=T(r,0)dr.A random variablek is used to decide whether to update the throughput time or not, which is either one or zero. The probability of k=1 equals wΔt and the probability of k=0 is 1-wΔt. Here, w is the update frequency. If k=1, the throughput time is updated and the new throughput time is given by a random number η(t), which is obtained from the throughput time distribution τ(r,WIP). If k=0, the current throughput time remains unchanged. In general, suppose that w depends on both η(t) and time t, then ((2.6b)) will be replaced by(2.6b′)P{k(t)=1}=Δtw(η(t),t),P{k(t)=0}=1-Δtw(η(t),t),dP{η(t)=r}=T(r,t)dr.
According to [18], w can be chosen as (2.7)w(r,t)=λ(t)rT-1,
where λ(t) is the influx of items and T-1=∫r-1T(r,t)dr. Based on the evolution (2.6a), (2.6b), and (2.6c), the algorithm procedure to solve the phase and TPT of an item for the discrete event and continuous production is then as follows.Step 1.
Set the initial phase to 0 as the item enters the factory, increase the WIPW(t) by one and adjust the probability distribution T(r,WIP) accordingly.Step 2.
Computew(r,t) based on (2.7), and sample the parameters k(t) and η(t) according to their distribution ((2.6b)).Step 3.
Compute the TPTτ(t+Δt) and the phase ξ(t+Δt) according to ((2.6b)).Step 4.
Update the distribution TPTT(r,WIP), according to the current WIP at the time t+Δt, and back to Step 2, repeat this loop until the phase ξ=1. Now the time when ξ=1 becomes the time of the item leaving the system.The functionξ(t) plays a role of a position which is advanced with a randomly changing velocity 1/η(t). It is, therefore, to derive a kinetic formulation. In the context of supply chains, these models are often referred to as traffic flow models [26]. A product arrives at the first stage at time t=a and moves from x=0 to x=1 along a trajectory with velocity 1/ηa(t), in such a way that it reaches the end x=1 at time t=e, then the velocity has to be chosen in the following way ∫ae(1/ηa(t))dt=1. The total WIP W(t) and the flux F(x,t) at any point x∈[0,1] are then obtained as follows: (2.8)W(t)=∫[H(ξa(t))-H(ξa(t)-1)]λ(a)da,F(x,t)=∫ηa(t)δ(x-ξa(t))λ(a)da,
where λ(a) is concentrated on the arrival times an for a discrete model or is a continuous function for a continuous production model.


## 2.2. Basic Partial Differential Equation (PDE) Models

Partial differential equation (PDE) models are actually continuum approximation of fluid models. Recently, PDE models for large-scale multiple reentrant production systems have become an important research topic. PDE models do have several advantages, that is, they might be preferred when evaluating the overall performance of large-scale manufacturing systems or solving the optimization problem for a longer time plan. Currently, the continuous models have been applied to many fields and have achieved some significant research results.This is appropriate to describe a semiconductor manufacturing fab involving a large number of items in many stages:ρ(x,t) is the density of the products with units (units/space) in the system, which is the conserved variable. Here, x denotes the degree of completion (DOC), x=0 describes raw products that have just entered into the factory, and x=1 denotes finished products that are ready to exit from the system. So, x is in the closed interval: x∈[0,1]. The total number of products in the system can be obtained by taking the integral of density ρ(x,t) of products over the stage variable x from 0 to 1. Therefore, the total WIP W(t) as a function of time can be obtained as follows: (2.9)W(t)=∫01ρ(x,t)dx.Assuming that there is a unique entry and exit for the system and the yield is 100%, PDE models can be given bellow according to the conservation law,(2.10)∂ρ(x,t)∂t+∂(v(ρ(x,t))ρ(x,t))∂x=0,x∈[0,1],t∈(0,∞),
where v(ρ(x,t)) is a velocity function that depends on the density ρ(x,t) only. For a reentrant supply chain system, we assume that v(ρ(x,t)) can be described by a state equation of the following form: (2.11)v(ρ(x,t))=v0(1-W(t)Wmax⁡),
where v0 is the velocity for the empty supply chain system, and Wmax is the maximal load (capacity of the supply chain system). Clearly, the velocity v(ρ(x,t)) is determined by the total WIP. The boundary condition for the start rate λ(t) of products entering the supply chain system at x=0 is then defined as(2.12)λ(t)=ρ(0,t)v(t).An arbitrary initial condition for the density of the products can be expressed as(2.13)ρ(x,0)=ρ0(x).The production process of the supply chain network can be described as an equivalentM/M/1 queue. The state variable ρeq denotes the equilibrium density of the supply chain system as a whole. Let ρ=ρeq, then ρ(0,t)=ρeq and v=veq=1/τ. Correspondingly, the associated cycle time in steady state is τ=1/veq. Since a job arriving at a queue with a processing rate is μ=1/vmax, the cycle time τ=μ(1+L) can be obtained according to the queuing theory, the equilibrium velocity therefore becomes(2.14)v(ρ)=vmax⁡1+∫01ρ(x,t)dx=vmax⁡1+w(t).Equation (2.14) is a widely used expression between v and ρ for large-scale multiple reentrant supply chain systems, notice that the velocity v only depends upon the WIP at stage x. Therefore, (2.10) can be rewritten as(2.15)∂ρ(x,t)∂t+v∂ρ(x,t)∂x=0,x∈[0,1],t∈(0,∞).Hence, assume that an initial WIP distributionρ0(x) in the factory is given, the resulting full PDE models for the single-product multiple reentrant supply chain systems are given as follows:(2.16)∂ρ(x,t)∂t+v(ρ)∂ρ(x,t)∂x=0,ρ(x,0)=ρ0(x),ρ(0,t)v(t)=λ(t),v(ρ)=vmax⁡1+∫01ρ(x,t)dx.If the influxλ(t) and the initial condition ρ0(x) are nonnegative, then the density will remain nonnegative. We use an upwinding scheme [27] to discretize the PDE which is given by the following equation:(2.17)ρ(xi,tj+1)=ρ(xi,tj)-ΔtΔxv(tj)[ρ(xi,tj)-ρ(xi-1,tj)],
where i=1,2,…,N, and j=0,1,…,M-1. Δt and Δx are the step sizes in time and space, respectively. Based on the boundary condition, the propagation scheme is given by(2.18)ρ(x0,tj+1)=ρ(x0,tj)-ΔtΔx[v(tj)ρ(x0,tj)-λ(tj)].Since the Courant-Friedrich-Levy (CFL) condition is necessary for stability, the time stepΔt and the space step Δx must satisfy the following formulation: |(Δt/Δx)vmax(t)|≤1. Here, vmax(t) is the maximum velocities in the system at time t.Based on the above formula, the density distribution of each momentρ(xi,tj) can be obtained. Therefore, the system throughput rate q(xN,tj) of each moment can also be computed as follows: (2.19)q(xN,tj)=ρ(xN,tj)v(ρ(xN,tj)).Furthermore, the total WIPw(t)=∫01ρ(x,t)dx can be obtained via the extended Simpson’s rule quadrature. Finally, the density distribution ρ(xi,tj) and throughput q(xN,tj) of each moment can be obtained.


## 3. The Modified Continuum Model

### 3.1. Numerical Experiment

Mini-Fab is a simplified model of the semiconductor production lines, which has all the important features of the reentrant semiconductor manufacturing systems, such as reentrant, different processing time, and batch production. Currently, many scholars have done a lot of research work based on the Mini-Fab. The Mini-Fab contains 5 machines grouped into 3 work centers, the product comprises of 6 processing steps and each work center has a reentrant step, which is shown in Figure2. The system designed in this study is FIFO (First In First Out).Figure 2
Process flow diagram of the Mini-Fab.For convenience, we make the following basic assumptions on the model.(1)
The product yield rate is 100%, namely, there is no rework problem.(2)
The system is a continuous production process for 24 hours a day.(3)
The system does not take into account the time of carrying, loading and discharging, adjusting equipment, premaintaining equipment, and downtime.Now assuming that there is one productD in the Mini-Fab, its processing steps and processing time are shown in Table 1.Table 1
Processing time of the productD at each step.
Machining centersProcessing time (hours)Machines A & BStep  1: 1.5Step  5: 1.5Machines C & DStep  2: 0.5Step  4: 1Machine EStep  3: 1Step  6: 0.5From the above table, we can see that the total processing time isP=0.25 (days). Based on the basic PDE model, the parameters are set as follows: vmax=4 (units/day), λ=5 (units/day). We typically start up with an empty production system, the total running time of the production lines is 10 (days). Let Δt=0.001 and Δx=0.01, then Δx and Δt that satisfy the CFL stability condition is formulated as follows: (Δt/Δx)vmax<1. The system throughput can be obtained through the simulation based on the basic PDE model of (2.16). Figure 3 indicates that the system throughput is about 3.2 (units/day) during the steady state.Figure 3
Throughput as a function of time for the PDE simulation.When the throughput of the basic continuous models is obtained, the corresponding Mini-Fab simulation model can be built using simulation package ExtendSim [28] to verify the validity of the models with a period of one year. The jobs are sorted by FIFO. From Figure 4, it can be seen that the throughput during the steady state is about 5 (units/day) through ExtendSim simulation.Figure 4
Throughput as a function of time for the ExtendSim simulation.Comparing Figure3 with Figure 4, obviously, it can be found that there are some differences between the throughput results of the two models. This is because the basic continuous model of (2.16) is built on the basis of a large number of materials and many reentrant steps in the systems and some important characteristics of semiconductor production systems are not captured. Therefore, the further investigation is needed to explore more precise models for multiple reentrant production systems.


### 3.2. The Modified Continuum Model

It is worth noting that the state equations could reflect the characteristics of systems—any changes of the multiple reentrant production systems may lead to a different state equation. Equation (2.14) is a quite general state equation, in order to capture some important features of multiple reentrant supply chain systems, a more specified relationship is required for computing the velocities numerically. Lefeber and Armbruster [29] presented a more sophisticated reentrant factory model through the use of integration kernels. For the general supply chain systems (nonreentrant systems), Sun and Dong [30] described several kinds of state equations and the corresponding cycle times. Although there are many kinds of continuum models currently, they do not reflect how the reentrant degree of the product impacts on the system performance. Hence, a new concept is introduced to describe the reentrant degree of production system.Definition 3.1.
Letα be a product reentrant factor, α equals to the ratio of the product’s processing time of reentrant steps and the product’s total processing time.Reentrant factorα is the property of product’s process flows, with the increase of reentrant factor α, the degree of reentrant of the product becomes larger and larger. Let P be the total processing time, P1 be the reentrant processing time, and P2 be the nonreentrant processing time. According to the above definition, the following formula can be obtained:(3.1)α=P1P1+P2=P1P.In reality, the velocity of products in the system is not only related with the WIP level, but also with the reentrant factor. Letw(t) be the WIP level in the system at a given time t, w(Δx,t) be the WIP level at interval Δx at time t, Δx be the interval of the completion of the product, and Δt be the processing time required to complete the interval Δx. Assuming that w(Δx,t) is proportional to Δt, and the system consists of two parts: reentrant processes and nonreentrant processes. Hence, αw(t) is the WIP level of the reentrant processes at a given time t, also referred as reentrant WIP. Similarly, (1-α)w(t) is the WIP level of nonreentrant processes at a given time t, also called as nonreentrant WIP.According to the queuing theory, the processing cycle time of the reentrant processτ1 can be expressed as follows:(3.2)τ1=[1+αw(t)]P1.As for the nonreentrant processes, assuming the total number of workstations ism, a1,a2,…,am are the corresponding processing times at each workstation, respectively, and τ2 is the nonreentrant processing cycle time, then we have(3.3)P2=a1+a2+⋯+am,τ2=∑i=1m(1+aiPw(t))ai=∑i=1mai+∑i=1mai2Pw(t)=P2+∑i=1mai2Pw(t).
So the total cycle time τ can be written in the following form: (3.4)τ=τ1+τ2=[1+αw(t)]P1+[P2+∑i=1mai2Pw(t)]=P+(α2+∑i=1mai2P2)w(t)P.
Then, the resulting new state equation for the velocity will be given by(3.5)v=1τ=1P+(α2+∑i=1mai2/P2)w(t)P.In order to avoid computational complexity of∑i=1mai2/P2, an approximate method is proposed to deal with the nonreentrant processes. Suppossing that WIP in nonreentrant processes follows the uniform distribution in each workstation, then the mean processing time P2/m and the mean WIP level at each workstation ((1-α)/m)w(t) can be obtained. According to the queuing theory, the mean cycle time of the nonreentrant processes at each workstation is [1+((1-α)/m)w(t)]P2/m, and the total cycle time of the nonreentrant processes τ2 is [1+((1-α)/m)w(t)]P2. Therefore,  the total cycle time τ can be expressed as follows:(3.6)τ=τ1+τ2=[1+αw(t)]P1+[1+(1-α)mw(t)]P2=P1+P2+[αP1+(1-α)mP2]w(t)=P+[α⋅αP+(1-α)⋅(1-α)Pm]w(t)=P+[α2+(1-α)2m]w(t)P.The corresponding new state equation can be obtained as follows:(3.7)v=1τ=1P+[α2+(1-α)2/m]w(t)P.Therefore, the resulting modified whole PDE model is given below:(3.8)∂ρ(x,t)∂t+v(ρ)∂ρ(x,t)∂x=0,ρ(x,0)=ρ0(x),ρ(0,t)v(t)=λ(t),v=1τ=vmax⁡1+[α2+(1-α)2/m]w(t).


### 3.3. Validity of the Modified Models

Once the new state equation is obtained, the validity of the modified models can be verified by the same example described in the previous section. According to Figure2, the total number of workstations is m=3,  and the steps  4 to 6 are reentrant steps. From Table 1, the reentrant factor of products α=0.5 can be easily obtained by definition. Assuming that the influx λ of items is 5 (units/day) and the total running time is 10 (days), the system starts running from an empty factory. With Δt=0.001, Δx=0.01, then Δx and Δt satisfy the CFL stability condition: vmax·Δt/Δx<1. Similarly, the modified continuous model (3.8) can also be solved via the upwind scheme, and the throughput of the modified continuum models can be obtained, which is plotted in Figure 5.Figure 5
Throughput as a function of time for the modified PDE model.It can be seen that the throughput is initially zero for the reentrant manufacturing systems. This is due to the time delay and the system initialization (the system starts up with an empty factory), then the system begins to have throughput at about 0.25 days and increases drastically until the throughput reaches a stable value of 5 (units/day) at about 1.5 days. Comparing Figure4 with Figure 5, it can be obviously found that the results are basically consistent with the simulation results obtained from ExtendSim. Therefore, it can be seen that the modified continuum models are effective for multiple reentrant supply chain systems.


## 3.1. Numerical Experiment

Mini-Fab is a simplified model of the semiconductor production lines, which has all the important features of the reentrant semiconductor manufacturing systems, such as reentrant, different processing time, and batch production. Currently, many scholars have done a lot of research work based on the Mini-Fab. The Mini-Fab contains 5 machines grouped into 3 work centers, the product comprises of 6 processing steps and each work center has a reentrant step, which is shown in Figure2. The system designed in this study is FIFO (First In First Out).Figure 2
Process flow diagram of the Mini-Fab.For convenience, we make the following basic assumptions on the model.(1)
The product yield rate is 100%, namely, there is no rework problem.(2)
The system is a continuous production process for 24 hours a day.(3)
The system does not take into account the time of carrying, loading and discharging, adjusting equipment, premaintaining equipment, and downtime.Now assuming that there is one productD in the Mini-Fab, its processing steps and processing time are shown in Table 1.Table 1
Processing time of the productD at each step.
Machining centersProcessing time (hours)Machines A & BStep  1: 1.5Step  5: 1.5Machines C & DStep  2: 0.5Step  4: 1Machine EStep  3: 1Step  6: 0.5From the above table, we can see that the total processing time isP=0.25 (days). Based on the basic PDE model, the parameters are set as follows: vmax=4 (units/day), λ=5 (units/day). We typically start up with an empty production system, the total running time of the production lines is 10 (days). Let Δt=0.001 and Δx=0.01, then Δx and Δt that satisfy the CFL stability condition is formulated as follows: (Δt/Δx)vmax<1. The system throughput can be obtained through the simulation based on the basic PDE model of (2.16). Figure 3 indicates that the system throughput is about 3.2 (units/day) during the steady state.Figure 3
Throughput as a function of time for the PDE simulation.When the throughput of the basic continuous models is obtained, the corresponding Mini-Fab simulation model can be built using simulation package ExtendSim [28] to verify the validity of the models with a period of one year. The jobs are sorted by FIFO. From Figure 4, it can be seen that the throughput during the steady state is about 5 (units/day) through ExtendSim simulation.Figure 4
Throughput as a function of time for the ExtendSim simulation.Comparing Figure3 with Figure 4, obviously, it can be found that there are some differences between the throughput results of the two models. This is because the basic continuous model of (2.16) is built on the basis of a large number of materials and many reentrant steps in the systems and some important characteristics of semiconductor production systems are not captured. Therefore, the further investigation is needed to explore more precise models for multiple reentrant production systems.


## 3.2. The Modified Continuum Model

It is worth noting that the state equations could reflect the characteristics of systems—any changes of the multiple reentrant production systems may lead to a different state equation. Equation (2.14) is a quite general state equation, in order to capture some important features of multiple reentrant supply chain systems, a more specified relationship is required for computing the velocities numerically. Lefeber and Armbruster [29] presented a more sophisticated reentrant factory model through the use of integration kernels. For the general supply chain systems (nonreentrant systems), Sun and Dong [30] described several kinds of state equations and the corresponding cycle times. Although there are many kinds of continuum models currently, they do not reflect how the reentrant degree of the product impacts on the system performance. Hence, a new concept is introduced to describe the reentrant degree of production system.Definition 3.1.
Letα be a product reentrant factor, α equals to the ratio of the product’s processing time of reentrant steps and the product’s total processing time.Reentrant factorα is the property of product’s process flows, with the increase of reentrant factor α, the degree of reentrant of the product becomes larger and larger. Let P be the total processing time, P1 be the reentrant processing time, and P2 be the nonreentrant processing time. According to the above definition, the following formula can be obtained:(3.1)α=P1P1+P2=P1P.In reality, the velocity of products in the system is not only related with the WIP level, but also with the reentrant factor. Letw(t) be the WIP level in the system at a given time t, w(Δx,t) be the WIP level at interval Δx at time t, Δx be the interval of the completion of the product, and Δt be the processing time required to complete the interval Δx. Assuming that w(Δx,t) is proportional to Δt, and the system consists of two parts: reentrant processes and nonreentrant processes. Hence, αw(t) is the WIP level of the reentrant processes at a given time t, also referred as reentrant WIP. Similarly, (1-α)w(t) is the WIP level of nonreentrant processes at a given time t, also called as nonreentrant WIP.According to the queuing theory, the processing cycle time of the reentrant processτ1 can be expressed as follows:(3.2)τ1=[1+αw(t)]P1.As for the nonreentrant processes, assuming the total number of workstations ism, a1,a2,…,am are the corresponding processing times at each workstation, respectively, and τ2 is the nonreentrant processing cycle time, then we have(3.3)P2=a1+a2+⋯+am,τ2=∑i=1m(1+aiPw(t))ai=∑i=1mai+∑i=1mai2Pw(t)=P2+∑i=1mai2Pw(t).
So the total cycle time τ can be written in the following form: (3.4)τ=τ1+τ2=[1+αw(t)]P1+[P2+∑i=1mai2Pw(t)]=P+(α2+∑i=1mai2P2)w(t)P.
Then, the resulting new state equation for the velocity will be given by(3.5)v=1τ=1P+(α2+∑i=1mai2/P2)w(t)P.In order to avoid computational complexity of∑i=1mai2/P2, an approximate method is proposed to deal with the nonreentrant processes. Suppossing that WIP in nonreentrant processes follows the uniform distribution in each workstation, then the mean processing time P2/m and the mean WIP level at each workstation ((1-α)/m)w(t) can be obtained. According to the queuing theory, the mean cycle time of the nonreentrant processes at each workstation is [1+((1-α)/m)w(t)]P2/m, and the total cycle time of the nonreentrant processes τ2 is [1+((1-α)/m)w(t)]P2. Therefore,  the total cycle time τ can be expressed as follows:(3.6)τ=τ1+τ2=[1+αw(t)]P1+[1+(1-α)mw(t)]P2=P1+P2+[αP1+(1-α)mP2]w(t)=P+[α⋅αP+(1-α)⋅(1-α)Pm]w(t)=P+[α2+(1-α)2m]w(t)P.The corresponding new state equation can be obtained as follows:(3.7)v=1τ=1P+[α2+(1-α)2/m]w(t)P.Therefore, the resulting modified whole PDE model is given below:(3.8)∂ρ(x,t)∂t+v(ρ)∂ρ(x,t)∂x=0,ρ(x,0)=ρ0(x),ρ(0,t)v(t)=λ(t),v=1τ=vmax⁡1+[α2+(1-α)2/m]w(t).


## 3.3. Validity of the Modified Models

Once the new state equation is obtained, the validity of the modified models can be verified by the same example described in the previous section. According to Figure2, the total number of workstations is m=3,  and the steps  4 to 6 are reentrant steps. From Table 1, the reentrant factor of products α=0.5 can be easily obtained by definition. Assuming that the influx λ of items is 5 (units/day) and the total running time is 10 (days), the system starts running from an empty factory. With Δt=0.001, Δx=0.01, then Δx and Δt satisfy the CFL stability condition: vmax·Δt/Δx<1. Similarly, the modified continuous model (3.8) can also be solved via the upwind scheme, and the throughput of the modified continuum models can be obtained, which is plotted in Figure 5.Figure 5
Throughput as a function of time for the modified PDE model.It can be seen that the throughput is initially zero for the reentrant manufacturing systems. This is due to the time delay and the system initialization (the system starts up with an empty factory), then the system begins to have throughput at about 0.25 days and increases drastically until the throughput reaches a stable value of 5 (units/day) at about 1.5 days. Comparing Figure4 with Figure 5, it can be obviously found that the results are basically consistent with the simulation results obtained from ExtendSim. Therefore, it can be seen that the modified continuum models are effective for multiple reentrant supply chain systems.


## 4. Numerical Experiments

Once the modified continuum models are obtained, in order to describe the multiple reentrant manufacturing systems from the micro- and macroperspectives, based on the DES and modified PDE models, a numerical experiment is carried out to demonstrate the consistency of the models for the Mini-Fab case. The WIP levels and outflux of the DES and the modified continuum models for the multiple reentrant supply chain systems are presented in this section.In this section, a mean field assumption is made first, namely, the TPT distributionT(r,t) is a given function while in reality it depends on the WIP W(t). The idea underlying this assumption is that, for many lots presenting in the system, the impact of one individual lot on the distribution T(r,t) is negligible. In engineering practices, the TPT distribution depends on the lot trajectories through the WIP W(t). Therefore, the lots can be treated as approximately statistically independent. The prescribed TPT distribution is shown as well as the upper and the lower limits as a function of WIP in Figure 6. The mean possible TPT is a constant between the maximal possible TPT and the lower possible TPT, while the minimally possible TPT remains constant independent of the WIP.Figure 6
Probability distribution of the TPT as a function of WIP.Same as the previous section, taking the multiple reentrant production system as an example, the system begins to run from an empty factory. The DES model is executed by generating a total of 200 lots. In order to facilitate the computation, we assume that the influx is defined as a constant value in the experiment first. The trajectories are computed by the random phase model, and the lots’ arrival timean is given by an+1=an+(1/λ(an)), where a0=0. The WIP levels and outflux of DES models are computed according to (2.8) with discrete arrival times an. For the modified PDE models, let vmax=4 (units/day), a uniform time step size Δt=0.001 and a uniform spatial interval length Δx=0.01, then Δx and Δt also satisfy the CFL stability condition: (Δt/Δx)vmax<1. Figures 7 and 8 show the WIP levels and outflux of the DES and modified PDE models, in which dotted line indicates the WIP levels and outflux variation of the DES models and the straight line represents the WIP levels and outflux variation of the modified PDE models, respectively.Figure 7
WIP levels of the DES and modified PDE models as a function of time.Figure 8
Outflux of the DES and the modified PDE models as a function of time.Comparing the WIP levels and outflux computed from the DES models to the solution of the modified PDE models, it can be observed that, although the WIP levels and outflux values of the two models are not exactly equal, the results are close enough to each other to substantiate the consistency of two models. The same results of the DES and modified PDE models are given for the steady state. The WIP levels of the two models show gradual increase, while the outflux of the two models first shows an increasing trend gradually to reach the maximum and then reaches a stable value.


## 5. Conclusions

In this paper, better analysis of the multiple reentrant manufacturing systems can be obtained if both micro- and macroperspectives are adopted. The discrete event simulation models are first proposed and their basic algorithm is also presented in detail. The low accuracy of the basic PDE models for the multiple reentrant supply chain networks is explained by a numerical experiment. In order to model such complex systems more precisely, a modified PDE model that takes into account the reentrant degree of the product is presented, while the validity of the modified PDE model is also illustrated through a numerical experiment for multiple reentrant supply chain systems. Then, based on the DES and modified PDE models, a numerical experiment is provided to compare the WIP levels and outflux changes. Meanwhile, some interesting observations are discussed. Once the results of the micro- and macro simulations are obtained, some analysis for the multiple reentrant manufacturing systems based on the multiscale methods becomes possible.



---

*Source: 139896-2011-05-03.xml*