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.

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. Figure 1 gives the structure of a reentrant manufacturing process. Here, 𝐡 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].

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: Section 2 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 of 𝑀 stages (𝑠1,𝑠2,…,𝑠𝑀), π‘Žπ‘šπ‘› denotes the time at which the lot number 𝑛 arrives at stage π‘ π‘š, and π‘’π‘šπ‘› denotes the time at which the lot number 𝑛 exits at stage π‘ π‘š (or the lot number 𝑛 arrivers at stage π‘ π‘š+1). πœπ‘šπ‘› denotes the processing time that the lot number 𝑛 takes at stage π‘ π‘š, which is known as the throughput time (TPT). The relationship between the arrival and exit times is given via the law(a)π‘’π‘šπ‘›=π‘Žπ‘šπ‘›+πœπ‘šπ‘›ξ€½πœ,(b)π‘‘π‘ƒπ‘šπ‘›ξ€Ύ=π‘Ÿ=π‘‡π‘šξ€·π‘Ÿ,π‘Žπ‘šπ‘›ξ€Έπ‘‘π‘Ÿ,(c)π‘Žπ‘›π‘š+1=π‘’π‘šπ‘›.(2.1) Here, π‘‡π‘š(π‘Ÿ,π‘Ž) denotes the time-dependent distribution of throughput times of stage π‘ π‘š, and 𝑃 denotes the probability distribution of the processing time. The TPT is determined by the state of the supplier π‘ π‘š at time π‘Žπ‘šπ‘› when lot number 𝑛 arrives. The throughput time distribution π‘‡π‘š is usually dependent on the total number lots π‘Šπ‘š(𝑑) processed at stage π‘ π‘š at time 𝑑, the so-called work in progress (WIP). The WIP π‘Šπ‘š(𝑑) is computed as π‘Šπ‘šξ“(𝑑)=π‘›π»ξ€·π‘‘βˆ’π‘Žπ‘šπ‘›ξ€Έξ€·βˆ’π»π‘‘βˆ’π‘Žπ‘›π‘š+1ξ€Έ,(2.2) where 𝐻 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 distribution π‘‡π‘š in (2.1) is written as π‘‡π‘š=π‘‡π‘š(π‘Ÿ,π‘Šπ‘š(π‘Žπ‘šπ‘›)). 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 π‘‡π‘š on the WIP π‘Šπ‘š 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 WIP π‘Šπ‘š(𝑑) and the influx πœ†π‘š(𝑑) 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 π‘Šπ‘š(𝑑) and the fluxes πœ†π‘š(𝑑) from stage π‘ π‘šβˆ’1 to π‘ π‘š. A random throughput time function 𝜏(𝑑) can be computed from the WIP obtained before. The WIP π‘Šπ‘š(𝑑) is then given by the conservation lawπ‘‘π‘Šπ‘‘π‘‘π‘š(𝑑)=πœ†π‘š(𝑑)βˆ’πœ‡π‘š(𝑑),πœ†π‘š+1(𝑑)=πœ‡π‘š(𝑑).(2.3) Here, πœ† and πœ‡ represent the influx and outflux of the lots, respectively. According to (2.1), the fluxes are given πœ‡π‘š(ξ€œπ‘‘)=𝛿(𝑠+𝜏(𝑠)βˆ’π‘‘)πœ†π‘š(𝑠)𝑑𝑠.(2.4)

For a large number of lots, the individual lots are replaced by a continuous quantity, and the fluxes πœ†π‘š are replaced by continuous functions. The WIP π‘Šπ‘š(𝑑) is evolved by the conservation law (𝑑/𝑑𝑑)π‘Š(𝑑)=πœ†(𝑑)βˆ’πœ‡(𝑑). Using the discrete time steps, the following formulation can be obtained: []π‘Š(𝑑+Δ𝑑)=π‘Š(𝑑)+Ξ”π‘‘πœ†(𝑑)βˆ’πœ‡(𝑑).(2.5)

So the question can be used to solve outflux πœ‡(𝑑) in terms of influx πœ†(𝑑).

However, if supplier π‘ π‘š is a reentrant system itself, then the lots arriving after time π‘Žπ‘šπ‘› will significantly influence the TPT πœπ‘šπ‘› since they will compete with lot number 𝑛 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 πœ‰(𝑑) and the throughput time 𝜏(𝑑) are introduced. The throughput time 𝜏(𝑑) is a random variable sampled from the distribution 𝑇(π‘Ÿ,𝑑). Given a prescribed TPT distribution function 𝑇(π‘Ÿ,𝑑), which is expressed in terms of WIP as 𝜏(π‘Ÿ,WIP)), the position and the throughput time after a small enough time interval Δ𝑑 are given byπœ‰(𝑑+Δ𝑑)=πœ‰(𝑑)+Ξ”π‘‘πœ‰πœ(𝑑),𝜏(𝑑+Δ𝑑)=(1βˆ’π‘˜)𝜏(𝑑)+π‘˜πœ‚(𝑑),(2.6a)𝑃{π‘˜(𝑑)=1}=𝑀Δ𝑑,𝑃{π‘˜(𝑑)=0}=1βˆ’π‘€Ξ”π‘‘,𝑑𝑃{πœ‚(𝑑)=π‘Ÿ}=𝑇(π‘Ÿ,𝑑)π‘‘π‘Ÿ,(2.6b)(0)=0,𝑑𝑃{πœ‚(0)=π‘Ÿ}=𝑇(π‘Ÿ,0)π‘‘π‘Ÿ.(2.6c)

A random variable π‘˜ is used to decide whether to update the throughput time or not, which is either one or zero. The probability of π‘˜=1 equals 𝑀Δ𝑑 and the probability of π‘˜=0 is 1βˆ’π‘€Ξ”π‘‘. Here, 𝑀 is the update frequency. If π‘˜=1, the throughput time is updated and the new throughput time is given by a random number πœ‚(𝑑), which is obtained from the throughput time distribution 𝜏(π‘Ÿ,WIP). If π‘˜=0, the current throughput time remains unchanged. In general, suppose that 𝑀 depends on both πœ‚(𝑑) and time 𝑑, then ((2.6b)) will be replaced by𝑃{π‘˜(𝑑)=1}=Δ𝑑𝑀(πœ‚(𝑑),𝑑),𝑃{π‘˜(𝑑)=0}=1βˆ’Ξ”π‘‘π‘€(πœ‚(𝑑),𝑑),𝑑𝑃{πœ‚(𝑑)=π‘Ÿ}=𝑇(π‘Ÿ,𝑑)π‘‘π‘Ÿ.(2.6bξ…ž) According to [18], 𝑀 can be chosen as 𝑀(π‘Ÿ,𝑑)=πœ†(𝑑)π‘Ÿπ‘‡βˆ’1,(2.7) where πœ†(𝑑) is the influx of items and π‘‡βˆ’1=βˆ«π‘Ÿβˆ’1𝑇(π‘Ÿ,𝑑)π‘‘π‘Ÿ. 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 WIP π‘Š(𝑑) by one and adjust the probability distribution 𝑇(π‘Ÿ,WIP) accordingly.

Step 2. Compute 𝑀(π‘Ÿ,𝑑) based on (2.7), and sample the parameters π‘˜(𝑑) and πœ‚(𝑑) according to their distribution ((2.6b)).

Step 3. Compute the TPT 𝜏(𝑑+Δ𝑑) and the phase πœ‰(𝑑+Δ𝑑) according to ((2.6b)).

Step 4. Update the distribution TPT 𝑇(π‘Ÿ,WIP), according to the current WIP at the time 𝑑+Δ𝑑, 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 πœ‰(𝑑) plays a role of a position which is advanced with a randomly changing velocity 1/πœ‚(𝑑). 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 𝑑=π‘Ž and moves from π‘₯=0 to π‘₯=1 along a trajectory with velocity 1/πœ‚π‘Ž(𝑑), in such a way that it reaches the end π‘₯=1 at time 𝑑=𝑒, then the velocity has to be chosen in the following way βˆ«π‘’π‘Ž(1/πœ‚π‘Ž(𝑑))𝑑𝑑=1. The total WIP π‘Š(𝑑) and the flux 𝐹(π‘₯,𝑑) at any point π‘₯∈[0,1] are then obtained as follows: ξ€œξ€Ίπ»ξ€·πœ‰π‘Š(𝑑)=π‘Ž(ξ€Έξ€·πœ‰π‘‘)βˆ’π»π‘Ž(ξ€œπœ‚π‘‘)βˆ’1ξ€Έξ€»πœ†(π‘Ž)π‘‘π‘Ž,𝐹(π‘₯,𝑑)=π‘Ž(𝑑)𝛿π‘₯βˆ’πœ‰π‘Ž(𝑑)πœ†(π‘Ž)π‘‘π‘Ž,(2.8) where πœ†(π‘Ž) is concentrated on the arrival times π‘Žπ‘› 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: 𝜌(π‘₯,𝑑) is the density of the products with units (units/space) in the system, which is the conserved variable. Here, π‘₯ denotes the degree of completion (DOC), π‘₯=0 describes raw products that have just entered into the factory, and π‘₯=1 denotes finished products that are ready to exit from the system. So, π‘₯ is in the closed interval: π‘₯∈[0,1]. The total number of products in the system can be obtained by taking the integral of density 𝜌(π‘₯,𝑑) of products over the stage variable π‘₯ from 0 to 1. Therefore, the total WIP π‘Š(𝑑) as a function of time can be obtained as follows: ξ€œπ‘Š(𝑑)=10𝜌(π‘₯,𝑑)𝑑π‘₯.(2.9)

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, πœ•πœŒ(π‘₯,𝑑)+πœ•π‘‘πœ•(𝑣(𝜌(π‘₯,𝑑))𝜌(π‘₯,𝑑))[]πœ•π‘₯=0,π‘₯∈0,1,π‘‘βˆˆ(0,∞),(2.10) where 𝑣(𝜌(π‘₯,𝑑)) is a velocity function that depends on the density 𝜌(π‘₯,𝑑) only. For a reentrant supply chain system, we assume that 𝑣(𝜌(π‘₯,𝑑)) can be described by a state equation of the following form: 𝑣(𝜌(π‘₯,𝑑))=𝑣0ξ‚΅1βˆ’π‘Š(𝑑)π‘Šmaxξ‚Ά,(2.11) where 𝑣0 is the velocity for the empty supply chain system, and π‘Šmax is the maximal load (capacity of the supply chain system). Clearly, the velocity 𝑣(𝜌(π‘₯,𝑑)) is determined by the total WIP. The boundary condition for the start rate πœ†(𝑑) of products entering the supply chain system at π‘₯=0 is then defined asπœ†(𝑑)=𝜌(0,𝑑)𝑣(𝑑).(2.12)

An arbitrary initial condition for the density of the products can be expressed as 𝜌(π‘₯,0)=𝜌0(π‘₯).(2.13)

The production process of the supply chain network can be described as an equivalent 𝑀/𝑀/1 queue. The state variable 𝜌eq denotes the equilibrium density of the supply chain system as a whole. Let 𝜌=𝜌eq, then 𝜌(0,𝑑)=𝜌eq and 𝑣=𝑣eq=1/𝜏. Correspondingly, the associated cycle time in steady state is 𝜏=1/𝑣eq. Since a job arriving at a queue with a processing rate is πœ‡=1/𝑣max, the cycle time 𝜏=πœ‡(1+𝐿) can be obtained according to the queuing theory, the equilibrium velocity therefore becomes𝑣𝑣(𝜌)=max∫1+10=π‘£πœŒ(π‘₯,𝑑)𝑑π‘₯max1+𝑀(𝑑).(2.14)

Equation (2.14) is a widely used expression between 𝑣 and 𝜌 for large-scale multiple reentrant supply chain systems, notice that the velocity 𝑣 only depends upon the WIP at stage π‘₯. Therefore, (2.10) can be rewritten asπœ•πœŒ(π‘₯,𝑑)πœ•π‘‘+π‘£πœ•πœŒ(π‘₯,𝑑)[]πœ•π‘₯=0,π‘₯∈0,1,π‘‘βˆˆ(0,∞).(2.15)

Hence, assume that an initial WIP distribution 𝜌0(π‘₯) in the factory is given, the resulting full PDE models for the single-product multiple reentrant supply chain systems are given as follows:πœ•πœŒ(π‘₯,𝑑)πœ•π‘‘+𝑣(𝜌)πœ•πœŒ(π‘₯,𝑑)πœ•π‘₯=0,𝜌(π‘₯,0)=𝜌0𝑣(π‘₯),𝜌(0,𝑑)𝑣(𝑑)=πœ†(𝑑),𝑣(𝜌)=max∫1+10.𝜌(π‘₯,𝑑)𝑑π‘₯(2.16)

If the influx πœ†(𝑑) and the initial condition 𝜌0(π‘₯) 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:πœŒξ€·π‘₯𝑖,𝑑𝑗+1ξ€Έξ€·π‘₯=πœŒπ‘–,π‘‘π‘—ξ€Έβˆ’Ξ”π‘‘π‘£ξ€·π‘‘Ξ”π‘₯π‘—πœŒξ€·π‘₯𝑖,𝑑𝑗π‘₯βˆ’πœŒπ‘–βˆ’1,𝑑𝑗,(2.17) where 𝑖=1,2,…,𝑁, and 𝑗=0,1,…,π‘€βˆ’1. Δ𝑑 and Ξ”π‘₯ are the step sizes in time and space, respectively. Based on the boundary condition, the propagation scheme is given byπœŒξ€·π‘₯0,𝑑𝑗+1ξ€Έξ€·π‘₯=𝜌0,π‘‘π‘—ξ€Έβˆ’Ξ”π‘‘ξ€Ίπ‘£ξ€·π‘‘Ξ”π‘₯π‘—ξ€ΈπœŒξ€·π‘₯0,π‘‘π‘—ξ€Έξ€·π‘‘βˆ’πœ†π‘—ξ€Έξ€».(2.18)

Since the Courant-Friedrich-Levy (CFL) condition is necessary for stability, the time step Δ𝑑 and the space step Ξ”π‘₯ must satisfy the following formulation: |(Δ𝑑/Ξ”π‘₯)𝑣max(𝑑)|≀1. Here, 𝑣max(𝑑) is the maximum velocities in the system at time 𝑑.

Based on the above formula, the density distribution of each moment 𝜌(π‘₯𝑖,𝑑𝑗) can be obtained. Therefore, the system throughput rate π‘ž(π‘₯𝑁,𝑑𝑗) of each moment can also be computed as follows: π‘žξ€·π‘₯𝑁,𝑑𝑗π‘₯=πœŒπ‘,π‘‘π‘—ξ€Έπ‘£ξ€·πœŒξ€·π‘₯𝑁,𝑑𝑗.(2.19)

Furthermore, the total WIP βˆ«π‘€(𝑑)=10𝜌(π‘₯,𝑑)𝑑π‘₯ can be obtained via the extended Simpson’s rule quadrature. Finally, the density distribution 𝜌(π‘₯𝑖,𝑑𝑗) and throughput π‘ž(π‘₯𝑁,𝑑𝑗) 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 Figure 2. The system designed in this study is FIFO (First In First Out).

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 product 𝐷 in the Mini-Fab, its processing steps and processing time are shown in Table 1.

From the above table, we can see that the total processing time is 𝑃=0.25 (days). Based on the basic PDE model, the parameters are set as follows: 𝑣max=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 Δ𝑑=0.001 and Ξ”π‘₯=0.01, then Ξ”π‘₯ and Δ𝑑 that satisfy the CFL stability condition is formulated as follows: (Δ𝑑/Ξ”π‘₯)𝑣max<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.

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.

Comparing Figure 3 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 𝑃 be the total processing time, 𝑃1 be the reentrant processing time, and 𝑃2 be the nonreentrant processing time. According to the above definition, the following formula can be obtained:𝑃𝛼=1𝑃1+𝑃2=𝑃1𝑃.(3.1)

In reality, the velocity of products in the system is not only related with the WIP level, but also with the reentrant factor. Let 𝑀(𝑑) be the WIP level in the system at a given time 𝑑, 𝑀(Ξ”π‘₯,𝑑) be the WIP level at interval Ξ”π‘₯ at time 𝑑, Ξ”π‘₯ be the interval of the completion of the product, and Δ𝑑 be the processing time required to complete the interval Ξ”π‘₯. Assuming that 𝑀(Ξ”π‘₯,𝑑) is proportional to Δ𝑑, and the system consists of two parts: reentrant processes and nonreentrant processes. Hence, 𝛼𝑀(𝑑) is the WIP level of the reentrant processes at a given time 𝑑, also referred as reentrant WIP. Similarly, (1βˆ’π›Ό)𝑀(𝑑) is the WIP level of nonreentrant processes at a given time 𝑑, also called as nonreentrant WIP.

According to the queuing theory, the processing cycle time of the reentrant process 𝜏1 can be expressed as follows:𝜏1=[]𝑃1+𝛼𝑀(𝑑)1.(3.2)

As for the nonreentrant processes, assuming the total number of workstations is π‘š, π‘Ž1,π‘Ž2,…,π‘Žπ‘š are the corresponding processing times at each workstation, respectively, and 𝜏2 is the nonreentrant processing cycle time, then we have𝑃2=π‘Ž1+π‘Ž2+β‹―+π‘Žπ‘š,𝜏2=ξ“π‘šπ‘–=1ξ‚€π‘Ž1+π‘–π‘ƒξ‚π‘Žπ‘€(𝑑)𝑖=ξ“π‘šπ‘–=1π‘Žπ‘–+ξ“π‘šπ‘–=1π‘Ž2𝑖𝑃𝑀(𝑑)=𝑃2+ξ“π‘šπ‘–=1π‘Ž2𝑖𝑃𝑀(𝑑).(3.3) So the total cycle time 𝜏 can be written in the following form: 𝜏=𝜏1+𝜏2=[]𝑃1+𝛼𝑀(𝑑)1+𝑃2+ξ“π‘šπ‘–=1π‘Ž2𝑖𝑃𝛼𝑀(𝑑)=𝑃+2+ξ“π‘šπ‘–=1π‘Ž2𝑖𝑃2ξƒͺ𝑀(𝑑)𝑃.(3.4) Then, the resulting new state equation for the velocity will be given by1𝑣=𝜏=1𝛼𝑃+2+βˆ‘π‘šπ‘–=1π‘Ž2𝑖/𝑃2𝑀(𝑑)𝑃.(3.5)

In order to avoid computational complexity of βˆ‘π‘šπ‘–=1π‘Ž2𝑖/𝑃2, 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 𝑃2/π‘š and the mean WIP level at each workstation ((1βˆ’π›Ό)/π‘š)𝑀(𝑑) can be obtained. According to the queuing theory, the mean cycle time of the nonreentrant processes at each workstation is [1+((1βˆ’π›Ό)/π‘š)𝑀(𝑑)]𝑃2/π‘š, and the total cycle time of the nonreentrant processes 𝜏2 is [1+((1βˆ’π›Ό)/π‘š)𝑀(𝑑)]𝑃2. Therefore, the total cycle time 𝜏 can be expressed as follows:𝜏=𝜏1+𝜏2=[]𝑃1+𝛼𝑀(𝑑)1+ξ‚Έ1+(1βˆ’π›Ό)π‘šξ‚Ήπ‘ƒπ‘€(𝑑)2=𝑃1+𝑃2+𝛼𝑃1+(1βˆ’π›Ό)π‘šπ‘ƒ2ξ‚Ήξ‚Έ(𝑀(𝑑)=𝑃+𝛼⋅𝛼𝑃+1βˆ’π›Ό)β‹…(1βˆ’π›Ό)π‘ƒπ‘šξ‚Ήξ‚Έπ›Όπ‘€(𝑑)=𝑃+2+(1βˆ’π›Ό)2π‘šξ‚Ήπ‘€(𝑑)𝑃.(3.6)

The corresponding new state equation can be obtained as follows:1𝑣=𝜏=1𝛼𝑃+2+(1βˆ’π›Ό)2ξ€»/π‘šπ‘€(𝑑)𝑃.(3.7)

Therefore, the resulting modified whole PDE model is given below:πœ•πœŒ(π‘₯,𝑑)πœ•π‘‘+𝑣(𝜌)πœ•πœŒ(π‘₯,𝑑)πœ•π‘₯=0,𝜌(π‘₯,0)=𝜌01(π‘₯),𝜌(0,𝑑)𝑣(𝑑)=πœ†(𝑑),𝑣=𝜏=𝑣max𝛼1+2+(1βˆ’π›Ό)2ξ€»./π‘šπ‘€(𝑑)(3.8)

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 Figure 2, the total number of workstations is π‘š=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 Δ𝑑=0.001, Ξ”π‘₯=0.01, then Ξ”π‘₯ and Δ𝑑 satisfy the CFL stability condition: 𝑣max⋅Δ𝑑/Ξ”π‘₯<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.

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 Figure 4 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 distribution 𝑇(π‘Ÿ,𝑑) is a given function while in reality it depends on the WIP π‘Š(𝑑). The idea underlying this assumption is that, for many lots presenting in the system, the impact of one individual lot on the distribution 𝑇(π‘Ÿ,𝑑) is negligible. In engineering practices, the TPT distribution depends on the lot trajectories through the WIP π‘Š(𝑑). 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.

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 time π‘Žπ‘› is given by π‘Žπ‘›+1=π‘Žπ‘›+(1/πœ†(π‘Žπ‘›)), where π‘Ž0=0. The WIP levels and outflux of DES models are computed according to (2.8) with discrete arrival times π‘Žπ‘›. For the modified PDE models, let 𝑣max=4 (units/day), a uniform time step size Δ𝑑=0.001 and a uniform spatial interval length Ξ”π‘₯=0.01, then Ξ”π‘₯ and Δ𝑑 also satisfy the CFL stability condition: (Δ𝑑/Ξ”π‘₯)𝑣max<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.

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.

Acknowledgments

The work presented in this paper has been supported by a Grant from the National High-Tech Research and Development Program (863 Program) of China (2008AA04Z104) and a Grant from National Natural Science Foundation of China (70871077).