Modeling and Analysis of Material Flows in Re-Entrant Supply Chain Networks Using Modiﬁed Partial Differential Equations

The basic partial di ﬀ erential equation (cid:3) PDE (cid:4) models for supply chain networks with re-entrant nodes and their macroscopic are proposed. However, through numerical examples, the basic continuum models do not perform well for multiple re-entrant systems. Then, a new state equation considering the re-entrant degree of the products is introduced to improve the e ﬀ ectiveness of the basic continuum models. The applicability of the modiﬁed continuum models for re-entrant supply chains is illustrated through numerical examples. Finally, based on the modiﬁed continuum model, numerical examples of di ﬀ erent re-entrant degrees are given, meanwhile, the changes in the WIP and outﬂux are analyzed in details for multiple re-entrant supply chain systems.


Introduction
In recent years, factories and production systems have become larger and more complicated. The re-entrant manufacturing system is a typical example, which indicates that work in process WIP repeatedly passes through the same workstation at different stages of the process routes. Figure 1 gives a simplified model of re-entrant systems. The model consists of three workstations and five buffers; the arrows in the figure indicate the processing path of a lot. From Figure 1, it can be seen that the lot passes through the same workstations M 1 and M 2 twice, which is a typical situation of a re-entrant system. A re-entrant supply chain network consists of many such manufacturing systems. Understanding the behavior of large supply chains under different policies and scenarios is a major issue for many businesses today. In large factories, no experiments can be done involving whole supply chains. Therefore, simulation models are developed, which substitute for the real environment. There are several methods to model the multiple re-entrant production flows: Petrinet, queuing network, fluid network, and partial differential equation PDE .
Petrinet is a mathematical tool to describe and analyze the logical relations of state changes for discrete event dynamic systems. It has been widely used in modeling, simulating, analyzing, and controlling of the discrete event dynamic systems. Compared with some other description tools, Petrinet model is especially easy to describe concurrent phenomena and simulate the parallel systems. On account of directly describing production processes, Petrinet model has been used to model the re-entrant manufacturing systems. However, it is extremely hard to solve the analytic solutions with the increase of system complexity. Wang and Wu 1 put forward an object-oriented hybrid Petrinet model of semiconductor manufacturing lines. Lin et al. 2 established a model of the re-entrant semiconductor production lines using Petrinets and studied on the stability of the system using bufferboundless approach. Dong and Chen 3 developed a modular modeling approach based on object-oriented predicate/transition nets OPTNs for the analysis of supply chain process models.
Queuing network model treats each workstation as an individual discrete queue, then these independent queues are connected to form queuing networks according to the production flow. Queuing network model can directly describe the process of production lines and supply chain networks. A methodology for supply chain inventory analysis and optimization was presented by linking production authorization PA strategy to queueing models 4 . The statistics of arrival flow of the fluid model in queuing systems has been studied 5, 6 . Dong and Chen 7 proposed a network of inventory-queue models for the performance modeling and analysis of an integrated supply chain network. S. Kumar and P. R. Kumar 8 focused on the analysis of queuing theory of the re-entrant manufacturing systems. However, the analytic solution could be obtained only for the small-size systems based on the queuing models. When the scale of production systems or supply chain networks become large, the number of system states would grow exponentially and a combinatorial explosion problem would occur. Therefore, the queuing models have some difficulties in analysis of the real large-scale re-entrant supply chain networks. Most of the queuing models can only be used to evaluate the stability of some scheduling policies and difficult to be used directly for the real-world supply chains.
Since most supply chains deal with individual parts and the processes that those parts undergo, the above two discrete event models can describe a lot of details of the studied system. While discrete event simulators have been highly successful to simulate single factories, they are computationally too expensive to simulate even a moderately complicated supply chains. Meanwhile, they are not scalable to a full supply chain. Therefore, a new modeling method called approximate modeling appears.
There are two kinds of approximate models: one is the so-called fluid model. It comes from traffic theory and was introduced by Newell 9 to approximately solve queuing problems. It treats the queue length l t as a continuous variable. The change rate of a queue length is equal to the arrival rate λ t minus the processing rate μ t when the queue length is not 0, otherwise the change rate of a queue length is 0. Hence, fluid models can quickly obtain the results by solving the ordinary differential equations. The time required for computing is not affected with the increase of material flows. The fluid models can be used not only to optimize the control system but also to reflect the long-term performance behavior of the system in the dynamic states.
However, one of the most important shortcomings of fluid models is that they cannot reflect the stochasticity of the production processes very well, that is, in the above equation, if λ t and μ t are mean rates, then this is a fully deterministic system, stochasticity is not modeled at all. Otherwise, if λ t and μ t are stochastic processes, then some theoretical analysis will be allowed. But this will degrade the advantages of a continuum models as a simulation tool 10 .
Dai and Weiss 11 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 12 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. The other approximate model is called partial differential equation PDE model. PDE models are actually continuum limit of fluid models. PDE models do have several advantages, that is, they are scalable; more detailed results can be found as compared to fluid models, and more importantly, they are amenable to optimization and control 13 .
Recently, the continuous models have been applied to many fields and have achieved some significant research results. Anderson 14 established the basic continuous model for supply chains and described production flow of the system using rate equations macroscopically. Lee et al. 15 studied the supply chain simulation with combined discretecontinuous modeling method. It integrates the wide applicability of the discrete event simulation DES and fast computation of the continuum models together. Armburster et al. 16 introduced the concept of materials' density and established the continuous model of large-scale re-entrant manufacturing systems and proposed new state equations so that the continuous model could be used to the real production systems. Compared with a DES model, van den Berg et al. 17 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. Unver et al. 18 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 to accurately describe the multiple re-entrant supply chain networks from macroscopic perspective, a modified model is proposed, which can reflect how the re-entrant degree of a product impacts the system performance. Based on the queuing theory, we restrict our attention on establishing a new state equation considering the re-entrant degrees of jobs or products, meanwhile, we will also verify the validity of the modified PDE model for multiple re-entrant semiconductor supply chain systems by 4 Journal of Applied Mathematics a numerical experiment. And finally, the changes of the WIP and outflux in the case of different re-entrant factors are analyzed based on the modified PDE model.
The structure of this paper is organized as follows: in Section 2, the basic continuum model is presented and a numerical experiment is carried out to verify the validity of the basic model. Section 3 introduces the concept of re-entrant factor to improve the basic continuum model and verifies its effectiveness through an example. In Section 4, based on the modified continuous model, the WIP profile and outflux are analyzed in the case of different re-entrant factors. Finally, some conclusions are given in Section 5.

Basic Model
Recently, continuum models for large-scale multiple re-entrant production systems have become an important research topic. Such a description is appropriate for a semiconductor manufacturing fab involving a large number of items in many stages. As ρ x, t is the conserved variable, it is the density of the products with units units/space in the system.
Here, x denotes the completion variable, 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, in the closed interval x ∈ 0, 1 . The total number of products in the system can be found 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: According to the conservation law, assuming that there is a unique entry and exit for the system and the yield is 100%, PDE models can be given by where v ρ x, t is a velocity function that depends on the density ρ x, t only. For a re-entrant supply chain system, we assume that v ρ x, t can be described by a state equation of the following form: Here, v 0 is the velocity for the empty supply chain system and W max 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 Journal of Applied Mathematics 5 An arbitrary initial condition for the density of the products can be expressed as The production process in the supply chain network is described as an equivalent M/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 v eq 1/τ. Correspondingly, the associated cycle time in steady state is τ 1/v eq . Since a job arriving at a queue with a processing rate is μ 1/v max , the cycle time τ μ 1 L can be obtained according to the queuing theory; therefore, the equilibrium velocity becomes v ρ v max Equation 2.6 is a widely used expression between v and ρ for large-scale multiple re-entrant supply chain systems; note that the velocity v only depends upon the WIP at stage x. Therefore, Equation 2.2 can be rewritten as Hence, assuming an initial WIP distribution ρ 0 x in the factory is given, the resulting full PDE model for the single-product multiple re-entrant supply chain systems is given as follows

2.8
If the influx λ t and the initial condition ρ 0 x are nonnegative, then the density will remain nonnegative. We use an upwind scheme 19 to discretize the PDE which is given by the following equation: 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 Since the Courant-Friedrich-Levy CFL condition is necessary for stability, the time step Δt and the space step Δx must satisfy the formulation as follows: | Δt/Δx v max t | ≤ 1.
Here, v max t is the maximum of all occurring velocities in the system at time t.
Based on the above formula, the density distribution of each moment ρ x i , t j can be obtained. Therefore, the system throughput rate q x N , t j of each moment can also be computed as follows: Furthermore, the total WIP w t 1 0 ρ x, t dx can be obtained via the extended Simpson's rule quadrature, then the following equation can be obtained:

2.12
Finally, the density distribution ρ x i , t j and throughput q x N , t j of each moment can be obtained.

Numerical Experiments
Mini-Fab is a simplified model of the semiconductor production line, which has all the important features of the re-entrant 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 tool sets; the product comprises of 6 processing steps, each tool set is visited twice, which is shown in Figure 2 20 .
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, and premaintaining equipment and downtime.

Machining centers Processing time hours Machines A & B
Step 1: 1.5 Step 5: 1.5 Machines C & D Step 2: 0.5 Step 4: 1 Machine E Step 3: 1 Step 6: 0.5 Now it is assumed that there is one product D in the Mini-Fab, processing steps and processing time shown in Table 1.
Based on the basic continuous model, the parameters are assumed as follows: P 0.25 days , v 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 Δt 0.001 and Δx 0.01, then Δx and Δt satisfy the CFL stability condition Δt/Δx v max < 1. The system throughput can be obtained through the simulation based on the basic PDE model of 2.8 . Figure 3 indicates that the system throughput is about 3.2 units/day during the steady state.
Once the throughput of the basic continuous models is obtained, the corresponding Mini-Fab simulation model can be built using simulation package ExtendSim 21 to verify the validity of the models with a period of one year. As can be seen from Figure 4 that the throughput during the steady state is about 5 units/day through the ExtendSim simulation.
Comparing Figure 3 with Figure 4, it is obviously found that there are some differences between the throughput results of the two models. This is because the basic continuous model of 2.8 is built on the basis of a large number of materials, and many re-entrant 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 re-entrant production systems.

Modified Model
It is worth noting that the state equations can reflect the characteristics of systemsany change of the multiple re-entrant production systems may lead to a different state equation. Equation 2.6 is a quite general state equation, in order to capture some important features of multiple re-entrant supply chain systems, so a more specified relationship is required for computing velocities numerically. Lefeber and Armbruster 22 presented a more sophisticated re-entrant factory model through the use of integration kernels. For the general supply chain systems non-re-entrant systems , Sun and Dong 23 described several kinds of state equations and the corresponding cycle times. Although there are many kinds of continuum models for describing the re-entrant production system currently, they do not reflect how the re-entrant degree of the product impacts the system performance. Hence, a new concept reflecting the re-entrant degree of the products is introduced.
Definition 3.1. Let α be a product re-entrant factor, α equals the ratio of the product processing time of re-entrant steps and the product total processing time.
Re-entrant factor α is the property of product process flows, with the increase of reentrant factor α, the degree of re-entrant of the product becomes larger and larger. Let P be the total processing time, let P 1 be the re-entrant processing time, and let P 2 be the nonre-entrant processing time. According to the above definition, we can obtain the following formula: In reality, the velocity of products in the system is not only related with the WIP level but also with the re-entrant factor. Let w t be the WIP level in the system at a given time t, let w Δx, t be the WIP level at interval Δx at time t, let Δx be the interval of the completion of the product, and let Δt be the processing time required to complete the interval Δx. Assume that w Δx, t is proportional to Δt and the system consists of two parts: re-entrant processes and non-re-entrant processes, hence, αw t is the WIP level of the re-entrant processes at a given time t and also referred to as re-entrant WIP. Similarly, 1 − α w t is the WIP level of non-re-entrant processes at a given time t and also called as non-re-entrant WIP.
According to the queuing theory, the processing cycle time of the re-entrant process τ 1 can be expressed as follows: As for the non-re-entrant processes, assuming the total number of workstations is m, a 1 , a 2 , . . . , a m are the corresponding processing times at each workstation, respectively, and τ 2 is the non-re-entrant processing cycle time, then we have 3.3 so the total cycle time τ can be written in the following form: Then, the resulting new state equation for the velocity will be given by

3.5
In order to avoid computational complexity of m i 1 a 2 i /P 2 , an approximate method is proposed to deal with the non-re-entrant process. Suppose that WIP in non-re-entrant processes follows the uniform distribution in each workstation, then the mean processing time P 2 /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 non-re-entrant processes at each workstation is 1 1 − α /m w t P 2 /m , and the total cycle time of the non-re-entrant processes τ 2 is 1 1 − α /m w t P 2 . Therefore, the total cycle time τ can be expressed as follows: 3.6 The corresponding new state equation can be obtained as follows:

3.7
Therefore, the resulting modified whole PDE model is given below: 3.8

Validity of the Modified Model
Once the new state equation is obtained, the same example in the previous section can be used to verify the validity of the modified models. According to Figure 2, the total number of workstations is m 3, the steps 4 to 6 are re-entrant steps. From Table 1, the products' re-entrant factor α 0.5 can be easily obtained by definition. Assume that the influx λ is 5 units/day and the total running time is 10 days , the system starts running from an empty state. With Δt 0.001, Δx 0.01, then Δx and Δt satisfy the CFL stability condition: v max ·Δt/Δx < 1. Similarly, the modified continuous model 3.8 can be solved via the upwind scheme, and the throughput of the modified models can be obtained, which is plotted in Figure 5.
It can be seen that the throughput is initially zero for the re-entrant system. 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 about 0.25 days and increases drastically until the throughput reaches a stable value of 5 units/day about 1.5 days. It is obvious that the results are basically consistent with the simulation results obtained from ExtendSim. Therefore, it can be seen that the modified models are more effective for multiple re-entrant supply chain systems.

Case Studies for the Different Re-Entrant Factors
In order to analyze the influence of the re-entrant degree on the WIP profile and outflux for the multiple re-entrant supply chain systems, several scenarios of the different re-entrant factors are presented in this section. Same as the previous section, let v max 4 units/day , a uniform time step size Δt 0.001 days , and a uniform spatial interval length Δx 0.01, then Δx and Δt also satisfy the CFL stability condition Δt/Δx v max < 1. We assume that the varying influx is prescribed first, as shown in Figure 6. Taking the multiple re-entrant production systems as an example, the system begins to run from an empty factory. We observe the changes of the WIP and outflux in the case of different re-entrant factors. As shown in Figure 7, black and red lines indicate the WIP and outflux, respectively.
From Figure 7, we can observe that, roughly, the WIP of the modified PDE model shows an increasing trend gradually first to reach the maximum and then it has a decreasing trend until it reaches the stable value. Meanwhile, the outflux of the modified PDE model rises slowly and then sharply reaches the maximum, after that the outflux shows a sudden decrease and then slowly reaches a stable value. Additionally, it is easy to see that the maximum WIP occures a small decrease first and then rises slowly with the increase of the re-entrant factor. Correspondingly, the maximum outflux appears a small increasing trend at first; when the WIP reaches a certain value, the maximum outflux becomes smaller and smaller with the increase of the re-entrant factor. These results are almost consistent with the actual situation; meanwhile, the materials flow can be better controlled according to the WIP profile and outflux.

Conclusions
In this paper, the basic continuum models for material flows are first proposed, and its low accuracy for the multiple re-entrant supply chain networks is explained by a numerical experiment. In order to model such systems more precisely, this paper presents a new state equation that takes into account the re-entrant degree of the product. The applicability of the modified continuum model for multiple re-entrant supply chain systems is also illustrated through a numerical example. Afterwards, based on the modified continuum models, some numerical examples on different re-entrant factors are provided, and the impacts of different re-entrant factors on the WIP profile and outflux changes are studied. Meanwhile, some interesting observations are discussed. The proposed model can be used to obtain more accurate results for material flows of the multiple re-entrant supply chain networks from macroscopic perspective.