Parameter extraction of complex production systems via a kinetic approach

Continuum models of re-entrant production systems are developed that treat 
the flow of products in analogy to traffic flow. Specifically, the dynamics of material flow through a re-entrant 
factory via a parabolic conservation law is modeled describing the product density and flux in the factory. The basic 
idea underlying the approach is to obtain transport coefficients for fluid dynamic models in a multi-scale setting 
simultaneously from Monte Carlo simulations and actual observations of the physical system, i.e. the factory. 
Since partial differential equation (PDE) conservation laws are successfully used for modeling the dynamical behavior of product flow in manufacturing systems, a re-entrant manufacturing system is modeled using a diffusive PDE. The specifics of the production process enter into the velocity and diffusion coefficients of the conservation law. The 
resulting nonlinear parabolic conservation law model allows fast and accurate simulations. With the traffic flow-like 
PDE model, the transient behavior of the discrete event simulation (DES) model according to the averaged influx, which 
is obtained out of discrete event experiments, is predicted. The work brings out an almost universally applicable tool 
to provide rough estimates of the behavior of complex production systems in non-equilibrium regimes.


1.
Introduction. In recent years, factories and production systems have become larger and more complicated. For this reason, a research endeavor has been initiated to find time and cost efficient ways of production for such supply chains. Moreover, a wide range of traffic flow theories and models have been developed to find the most effective means of production. While discrete event simulators have been highly successful in simulating single factories, they are computationally too expensive to simulate even a moderately complicated production system [1]. Our goal in this paper is to present a computationally efficient way to simulate stochastic production systems.
A supply chain is a network of organizations comprising suppliers, manufacturers, customers, in some cases even distributors and retailers to dispatch any type of goods, services, and information to customers enabling the members gain competitive advantage in the market place. To gain competitive advantage, a supply chain should function efficiently and effectively to maximize the profits of the members and minimize the delivery time of goods. The performance of the whole chain depends on how efficient and effective each member functions and in the literature there has been a host of research on efficient and effective functioning of supply chains. Higher customer satisfaction and profits, significant determinants of the performance of supply chains, might be obtained by decreasing conflicts within the chain and creating superior customer value.
Understanding the behavior of large supply chains under different polices 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. Especially in recent years, fast scalable simulations of production flows in a supply chain have become a very important research topic. The long term goal of a supply chain simulation is to optimize and control the production across the whole supply chain. Since most production deals with individual parts and the processes that these parts undergo, discrete event simulators would be the regular method of choice for accurate simulations.
While discrete event simulators have been highly successful in simulating single factories, they are computationally too expensive to simulate even a moderately complicated production system. Also, they are not scalable to a full supply chain [1]. Alternative models that endow supply chain nodes with fixed production capacitated systems and fixed lead times are not accurate enough since they do not take into account the fact that capacitated system responds nonlinearly to increases in demand close to the limit of the production capacity.
The approach followed in this paper is to extract throughput times (TPTs) and work in progress (WIP) levels from a discrete event simulation, replacing actual observations of the physical system for the purposes of this work, and then to transform them into transport coefficients of a macroscopic conservation law. Using discrete event simulations we estimate the transport coefficients; namely velocity (convection, C) and diffusion (D). Velocity coefficient stands for the local velocity of each particle throughout the system. Diffusion coefficient exists due to the reentrance of the particles in the system, which will be similar to going backwards in the stages. Since we deal with stochastic systems; i.e., sometimes some machines might be down, this will cause both velocity and diffusion coefficients to be varying throughout the system. Ultimately, after extracting C and D, we are able to solve the diffusive PDE numerically. The PDE model that we use in this project is the conservation law given as follows: where x and t are space and time variables respectively. In the computations, x is considered as "stage" of the processes, which is actually an artificial "space" (road) in the simulations [13].
Here, ρ(x, t), which is actually the local work in progress, denotes the density of the particles and F (x, t) denotes the flux, both of which are used to model the dynamics of material flow through the re-entrant factory. Density of products (ρ(x, t)) is the conserved variable in the system with units [parts/space]. The flux F (x, t) is the number of parts per time with units [parts/time]. The velocity and diffusion coefficients are extracted from observations of the system (replaced by a discrete event simulation model for the purposes of this paper). Since PDE models are amenable to optimization and control [10], we basically merge the randomness and optimization by using both DES and PDE models.
We base our theory on the assumption that we have the following type of data available. The production process consists of M stages. a m n denotes the time at which lot number n arrives at stage m. Thus, the basic input of our macroscopic model consists of a table of arrival times a m n , n = 1 : N , m = 1 : M . From this table, the throughput times τ m n for each stage as well as WIP in any part of the supply chain can be easily computed. For instance, given the arrival times a m n the WIP in the m-th supplier W m (t) is computed as where H denotes the usual Heaviside function, which is also known as the integral of the Dirac delta function. As lots/particles go through the stages, we compute their velocities and variances depending on observed parameters. Since each particle has a different processing time at each stage, it has a varying velocity through out the system, so this causes a variance in velocity for each particle. According to fluid dynamics theory, variance results in diffusion. This results in computing TPTs, means, and variances depending on a chosen set of state parameters. We record the arrival times when lots pass through each stage in the production process. At the same time we record some macroscopic quantity (such as the WIP) of the system at each of the arrival times a m n , which gives a statistical distribution of throughput times, parameterized by this macroscopic quantity. We extract the transport coefficients; namely, velocity (convection) and diffusion coefficients of the particles in the production system, using DES. Detailed information on all features of DES can be found in [5] and [8]. The velocity and diffusion coefficients in (1) are then related to the mean and variance of this throughput distribution. Here, our goal is to form a compact model which is amenable itself to optimization. PDE models have the important advantage of being fast and that optimization can be applied to obtain accurate results. Our ultimate goal is to predict and optimize the transient behavior of a given system from observation without knowing details about the internal mechanisms of the system; i.e. we always want to know how the system behaves even though the features of the system is changed; for example, if a processor is changed/canceled in the system or if a machine is down for a while...etc. Our ultimate goal is to predict and optimize the transient behavior of a given system from observation without knowing details about the internal mechanisms of the system; i.e. we always want to know how the system behaves even though the features of the system is changed.
The contents of the rest of the paper are as follows: In Section 2, we present the theoretical background of traffic flow-like PDE models and how to relate TPTs and WIPs to velocity and diffusion coefficients. In Section 3, we explain the model hierarchies that we used to derive the flux function and how we compute the TPT distribution and its mean and variance from the given data. Finally, in Section 5, we apply the theories that we cover in this paper on a model problem and give numerical results comparing the discrete event experiments to the PDE model.

2.
Traffic flow-like PDE models. Traffic flow-like PDE models are continuum limits of fluid models. These models have several advantages: They are scalable and, most important of all, they are amenable to control and optimization [10]. There are two different types of PDE modeling: First principle modeling and direct modeling from observations. In the first principle modeling; on average there exists a functional relation between the flux and the WIP, which is known as a clearing function. Whereas in direct modeling from observation, which we do in our PDE modeling, the clearing function is replaced by a probability distribution function that is obtained from the observed data (DES in our case). After we obtain the data (randomly generated times) from DES, we get the TPT's and then, we extract the transport coefficients by using mean and variance of TPT's. We can then solve the PDE-conservation law given in Equation (1).
In our traffic flow-like PDE model, we introduce an artificial variable that represents each product's percentage of completion, or stage. The stage x goes from 0 (0%) to 1 (100%). x = 0 denotes the start of a product into the factory, in other words, the product enters the factory, and x = 1 denotes the end of a product that means the product leaves the factory. The arrival rates of the products at stages yield a flux function. Both the density and flux functions are very important in PDE modeling such that we use them to model the dynamics of material flow through the whole factory.
The total number of individual products in the system at time t that are at most 100x% of the way through the sequence of stages (machines) that those products must visit before leaving the system can be found at each time t ≥ 0 as follows The PDE model that we use in this project is a parabolic conservation law. A conservation law is a relation asserting that a specific quantity is conserved. Basically, for a quantity, to be conserved means that whatever enters to the system has to come out of that system after some time. For instance, in a manufacturing system, all the products that enter the system will exit the system. So, the quantity of the products will be conserved.
In our setting, a conservation law can be defined as a partial differential equation that expresses the fact that some physical quantity is locally conserved in a fluid or other continuous physical system, such as energy, momentum, or the quantity of fluid itself. Conservation law in the integral form is given by The integral in Equation (3) gives the number of products per time between x = a and x = b.
3. Hierarchy of models. We present a model hierarchy to derive the flux function F given in Equation (1) for a given production system. The model hierarchy consists of three parts, namely; a particle model for the trajectories of individual parts, a kinetic model for the phase space density, and a diffusion convection equation.
Particle and kinetic models were explained in detail in [3], however brief description of these models are also given in this paper below.
3.1. Particle models. In particle models, the evolution of a large ensemble of lots (particles) is modeled by describing the trajectories of each individual lot in phase space by a set of Newton equations. For a continuous time model, we introduce the trajectories as x = ξ(t) and v = ∂ξ ∂t . The velocity v is updated periodically from a velocity distribution V, which we get through the TPT's obtained from DES, so that we obtain a random variable v from where P{ 1 τ m n = v} denotes the probability of a local velocity 1 τ m n to be equal to v and V m denotes the velocity distribution at m th stage. Here, n and m stand for the numbers of lots and stages respectively; and a m n represents the time for each of these lots at each stage. Now, the position ξ(t) after a small time interval dt is given by: where v is the random velocity variable. Similarly, we have a rule for updating the velocity v according to a random variable κ, which we calculate by flipping a coin. The probability of κ = 1 is wdt and the probability of κ = 0 is 1 − wdt. Here w is the update frequency. In Section 3.2, we explain how to choose the update frequency w. Now, according to this rule, if κ = 0, we maintain the current velocity, and if κ = 1, we update the velocity v from the probability distribution (4) given by the formula Here, v is also a random velocity variable that we get from the velocity distribution V, which can be seen in Equation (4). Now, we describe in detail how to compute the velocity distribution V that is obtained from discrete event experiments. As we run a simulation, or observe an actual system, where N parts pass through M stages s 1 , ..., s M , we record the times a m n when part number n arrives at stage s m and finally at the exit, which corresponds to stage s M +1 . At the same time, we record some macroscopic state variables, denoted by Z m n , at every time a m n . Z m n = (Z 1 (n, m), .., Z K (n, m)) is in general a vector of K components. These Z components can be variables like the total WIP, downstream WIP (at higher stages), the upstream WIP (at lower stages), or just the station index m itself. We eventually form two tables: where there are different possibilities that the macroscopic state variable Z might depend on. Z m n = m if we want to make the velocity dependent on space only, or Z m n = W (a m n ) for dependence on the total WIP only, or Z m n = (W <m (a m n ), W >m (a m n )) for dependence on the WIP in front and behind the part. From one of these possibilities, we compute a distribution for the velocities, dependent on the global variables Z. We map the work stations , the collected data give a discrete probability distribution dependent on the macroscopic state variable Z of the form For the final model used in this work, the transport coefficients needed will be given solely in terms of these means and variances [12].
We are dealing with highly re-entrant production systems with many machines and buffers. Throughout the factory, machines with different speeds process the products. This accounts for the fact that in a re-entrant system, later arrivals influence the processing time of a particular product.
Equations (5) and (6) for the position and velocity define a stochastic process for the evolution of the products. So, briefly, the update algorithm for position (ξ) and velocity (v) is given by (c) ξ(0) = 0 As it can be seen in Equation (8), for every ∆t, we choose a random number κ, which is either zero or one, to decide whether to update the velocity v or not. We update the velocity at every time step of length ∆t with probability ω∆t. In general we let ω depend on v(t) and ξ(t). So, (8)(b) will be replaced by Choosing update frequency ω is discusses in Section 3.2. For analytic purposes, we would prefer to have a continuous time model, where the function ξ(t) is not given only at discrete time intervals. This could be achieved by the limit ω → ∞, giving a stochastic process. However, in this limit we would lose the information about the update frequency ω, which is a key parameter for controlling the stochasticity of the model.

3.2.
Kinetic models. This section is leading to a kinetic equation for the evolution of the density of the products at stage x with velocity v.
The update algorithm (8) is similar to a Monte Carlo method for the solution of a Boltzmann equation (see [7] for details). The Boltzmann equation is the kinetic equation for the time (t) evolution of the kinetic density function f (x, v, t) in oneparticle phase space, where x and v are position and velocity, respectively.
The function ξ(t), defined in Section 3.1, plays the role of a position which is advanced with a randomly changing velocity v(t). It is therefore tempting to derive a kinetic formulation. These models are often referred to as traffic flow models [2], [9], because of the analogy to lots moving on a freeway from start to completion. A product arrives at the first stage at time t = a and moves from x = 0 to x = 1 along a trajectory with velocity v 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 such a way that e a v a (t) dt = 1 holds.
The particle model for a product that arrives at the first stage of the system at t = a together with the phase speed model (8) is given as follows: The trajectory ξ a (t) satisfies the same Newton equations as in (8) together with the initial conditions The total WIP W (t), given in terms of the WIP density ρ(x, t), and the flux F (x, t) at any point x ∈ [0, 1] are then given by Here λ(t), the influx density, is concentrated on the arrival times a 0 n of parts in the system, i.e.
holds. From the definition given in (10), it is obvious that the WIP density ρ and the flux F will satisfy the conservation law We make a mean field assumption for the rest of this section, that is we will assume that the velocity distribution V(v, x) is a given function of velocity and stage. In reality, the velocity distribution will depend on the lot trajectories ξ through the WIP W , given by (10). The idea underlying this assumption is that, for many lots present in the system, the impact of one individual lot on the distribution V is negligible. Therefore, the lots can be treated as approximately statistically independent. This is a standard approach for a system of many molecules in gas dynamics [7]. Making this mean field assumption, we have the following theorem: Let f a (x, v, t) denote the probability density for finding ξ a , v a , defined by (9), at position x and velocity v, i.e. dP{ξ a (t) = x, v a (t) = v} = f a (x, v, t) dxdv. Then, f a satisfies the initial value problem in the limit ∆t → 0, weakly in x, v and t.
Using Theorem 3.1, we can derive the corresponding initial boundary value problem for the kinetic density f (x, v, t) of the number of lots at position x with velocity v.
Proof. The proof consists of essentially summing up all possibilities for the values of κ(t), ξ a (t) and v a (t) and weighting them with their respective probabilities. This takes the form Integrating with respect to x , v and u in the above two integrals gives where we have used the fact that V(u, x) du = 1 holds. We now multiply (13) by a test function ψ(x, v, t) and integrate over t > a giving where H denotes the Heaviside function. Shifting the arguments in the integrals gives The left hand side of the above equation can be rewritten as (15) Expanding the right hand side of (14) gives Setting this equal to (15) and dividing by ∆t gives which for ∆t → 0 is the weak formulation of (12). Now, we derive the corresponding initial boundary value problem for the kinetic density f (x, v, t). We define f by (16) f is not a probability density since the trajectory ξ a will only exist for t > a and the number of lots will not be constant over time. The lot density f will satisfy the same Boltzmann equation as the probability density f a with a more natural influx boundary condition replacing the initial condition (12)(b). Now, we have the following theorem: The lot density f , given in (16), satisfies the initial boundary value problem (IBVP) weakly in x, v and t. The Boltzmann equation (17) is a generalization of the kinetic model derived in [2]. The random phase approach has introduced the concept of collisions, i.e. the integral operator on the right hand side of (17), into the particle model.
Proof. The weak form of the initial value problem (12) implies that for any test function ψ(x, v, t) the relation Now, f a (x, v, t) and λ(a) vanish identically for x < 0 and a < 0, respectively. Therefore, f (x, v, t) vanishes identically for x < 0, t < 0. Thus, (19) is the same as (20) is the weak formulation of the boundary value problem (17).

3.2.1.
Choosing the update frequency. Now, it remains to choose the update frequency ω(v, t). We remark that there necessarily is some arbitrariness in this choice.
The discrete event simulation model yields one particular realization of a stochastic process. To compute meaningful answers, one therefore will run a series of discrete event experiments and compute averages. The more experiments, the closer the result will be to the average, where a degenerate throughput time distribution (a δ function concentrated at the mean) is used. So, the size of ω has to be chosen in the same way as the number of realizations of the discrete event simulation model is chosen. The guiding principle here has to be that the overall number of random numbers generated is the same. The basic idea of the approach presented so far is to modify the discrete event simulation model in order to take into account the transient changes in throughput times and influx, which will influence the outflux in the case of re-entrant nodes. So, on the other hand, if there are no changes in influx and throughput time, the discrete event model should correspond to a solution of the Boltzmann equation (17). Translating the discrete event simulation into kinetic model along the lines of the approach in [2] gives the following picture. For an arrival time a, we pick a velocity v according to the distribution V(v, x) and move with the constant velocity v = 1 τ until reaching x = 1 at time t = a + τ . This gives the relations for the lot trajectory in phase space by We denote by f d a (x, v, t) the corresponding probability density in phase space, i.e.
. The probability density f d (x, v, t) corresponding to (16) is then given by For the particle model to be equivalent to the discrete event simulation model, f d has to be a solution of the IBVP (17). Since f d is a function of t − x v , it satisfies ∂ t f d + v∂ x f d = 0. f d also satisfies the boundary condition (17)(b). Therefore, the residual obtained from inserting f d into (17) is given by the collision term and Q[f d ] = 0 has to hold in the case that neither V nor λ depend on time, giving which implies ω(v, t) = vγ(t) for some arbitrary function γ(t). This makes sense since we update the phase speed more frequently for shorter current throughput times, i.e. for a faster transient response of the system. So, with this choice, the discrete event simulation and the random phase model (8) are equivalent for constant V and λ and any function γ(t). In the general case, we choose γ(t) such that we roughly compute the same number of random numbers for either model. Using the discrete event simulation, we pick ∆tλ(t) random numbers for the velocity in the time interval ∆t. Using the random phase model, we pick ∆tγ(t) vV(v, x)dv in the same time interval. w and V 1 (the first velocity component from the distribution V) are now functions of v and x. Therefore, we set the update frequency ω(v, x) to where λ 0 is a characteristic value of the influx density λ(t).

Diffusion convection equation.
In this section, we explain how we obtain the velocity (convection) and diffusion coefficients. We express the flux F in terms of the density ρ via a functional expansion, namely the Chapman -Enskog expansion [6]. The Chapman -Enskog expansion is an essential asymptotic tool of the Boltzmann equation (17) for a regime where average throughput time and the time intervals between velocity updates in the particle model are small compared to the over all time scale. Fluid equations for the simulation of throughput times have been studied in [2] and [4], leading to simple convection equations for the density of lots as a function of degree of completion (DOC). We essentially recover these results but obtain as a generalization, an additional diffusive term arising from the variance of the velocity distribution V. A fluid dynamic approximation to the Boltzmann equation (17) is obtained by taking integrals with respect to the velocity v. Defining the density ρ(x, t) and the flux F (x, t) by we obtain the conservation law Here, ρ and F , defined by (22), represent expectation values of the position and flux densities, respectively. The goal is now to find an expression for the flux F in terms of the density ρ using an asymptotic solution of the Boltzmann equation (17) large time scales. We start by scaling the involved probability distributions and by bringing the Boltzmann equation (17) into an appropriate dimensionless form.
From (17), we have that the IVBP for the unscaled Boltzmann equation for the part density is of the form The appropriate scale for the velocity (and therefore also the time scale) are given by the probability distribution V(v, x). We choose the velocity scale v 0 = i.e. as the spatial average of the mean of the velocities, and rescale the probability distribution V correspondingly: (We denote the scaled and dimensionless variables with subscript s here.) We choose the over all time scale to be 1 v0 and define the scaled time as t = ts v0 . We scale the influx density λ(t) by its characteristic value λ 0 and scale the scattering frequency ω accordingly, setting making λ s and ω s dimensionless O(1) quantities. Finally, we scale the kinetic The scaled initial boundary value problem for the Boltzmann equation (17) now reads: The diffusion convection equation will be derived in the limit ε → 0. Note that, according to Section 3.1 the average expectation v 0 is of order O( 1 M τ ). Thus, the timescale t 0 = 1 v0 = O(M τ ) corresponds to the total average throughput time through the whole chain, and ε = v0 λ0 << 1 implies that the intervals between arrivals in the chain are small compared to the whole TPT, i.e. there are many parts in the whole system. Now, using the Chapman -Enskog expansion, the flux F is expressed for ε << 1 in terms of the density ρ [6]. For notational simplicity, we drop the subscript s from here on.
In the limit ε → 0, equation (24) will be dominated by the right hand side. We therefore define the collision operator Q[f ] as and write the Boltzmann equation (24) as The relation between the density ρ and the flux F in (22) is now obtained via asymptotics for ε → 0, in other words by a Chapman -Enskog procedure. In general, the Chapman -Enskog expansion consists of splitting the dynamics of the Boltzmann equation (26) into a slow evolution on the kernel of the operator Q and a fast evolution on its orthogonal complement. While the Chapman -Enskog expansion is in general a nonlinear procedure [6], its complexity reduces significantly in the case of linear collision operators, and it can be written in terms of projection operators. To this end, we first define the kernel of the operator Q in (25). The definition of Q is the following: When we apply Q on a function, we take the integral of that function and multiply by probability distribution V(v), and then subtract that function from the multiplication obtained. The operator Q given in (25) is a Bhatnagar-Gross-Krook (BGK) type collision operator, because the probability of the final state (v f ) is independent of the probability of the initial state (v i ). BGK is a simple Computational fluid dynamics technique used in the Boltzmann equation. This simple form of the collision operator allows for a relatively simple analysis of the kernel of Q.
Obviously, kernel elements are given by multiples of the function V(v,x) ω(v,x) . Given the form (24)(c) of the scattering frequency ω, we define the normalized kernel element as A projection on the kernel is then given by the operator P is a projection operator since the kernel density W(v, x) is normalized, i.e. W(v, x) dv = 1 ∀x ∈ W, and therefore P 2 = P holds. For the following, it will be notationally convenient to define symbols for the moments of the distributions V and W. We define Because of (27), the moments V j and W j are obviously related. They satisfy the relations We have the following theorem.
Theorem 3.3. For the dimensionless parameter ε << 1, the flux F in (23) is given in terms of the density ρ as with the velocity (convection) coefficient C and the diffusion coefficient D given by where the moments W j of the kernel density W are defined as in (29). The diffusion coefficient D(x) satisfies D ≥ 0.
Proof. We split the function space B for the kinetic density function f into the linear hull of the kernel and its complement: Comparing coefficients gives The diffusion coefficient D is always non -negative, since it can be written as a variance, i.e. we have Hence, we obtain the convection and diffusion coefficients and the asymptotic conservation law given by (36) The conservation law (36) can be used for modeling complex re-entrant production systems successfully.
4. Parameter extraction. In this section, we describe how to compute the transport coefficients, i.e. the convection and the diffusion coefficients from the velocity distribution V (see Section 3.1 for the derivation) for the macroscopic model -dependent on stage x and time t.
According to the fluid dynamics theory, the mean and the variance of the distribution V depending on the macroscopic state variable Z, obtained from the observations (DES), are related respectively to the average velocity coefficient C(x) and the diffusion coefficient D(x) in a fluid dynamics approximation of the form where ρ = ρ(x, t) represents density of the products dependent on stage and time, and F (ρ(x, t)) represents the product flow (flux).
Using the theory developed in Section 3.3, we relate the velocity (C) and the diffusion (D) coefficients to the mean and variance of the distribution V. After normalization, formulas for the velocity (C) and the diffusion (D) are respectively given by where 1 M denotes the length of space interval from one stage to the other and Ω = v 2 m σ 2 m denotes the variation coefficient, which is a measure of stochasticity. As Ω gets large, the stochasticity increases; and as Ω gets smaller, the system becomes more deterministic. The derivation of velocity coefficient is obviously obtained from the division of average distance by average time; whereas diffusion coefficient is derived from the particle model by long time averaging, which is known as "Chapman -Enskog expansion". The proof and the details of the derivations of the formulas can be found in Section 3.3.
Note that the PDE given in (1) represents in general a nonlinear problem since the transport coefficients C and D depend on the state variable Z which, in turn, will depend on the density ρ.

5.
A model problem. We carry out a numerical experiment in order to compare a complex system to the PDE model. The purpose of this comparison is to substantiate the agreement in terms of WIP levels and fluxes as functions of time at each machine stage.
As a semi -realistic example, we use a simplified model of a semiconductor fab. This model consists of nine machine clusters, each performing tasks such as etching, lithography and deposition. Each cluster consists of multiple machines. So, there is a total of 200 machines in the system. The simplified production process consists of 26 stages, i.e. each part has to visit an individual cluster multiple times, and we are dealing with a highly re-entrant production system. The re-entrant structure in semiconductor manufacturing arises from the fact that chips are produced in layers. So, each time a layer is deposited, the same wafer has to go through the sequence of lithography -deposition -etching again.
For a detailed explanation of a simplified semiconductor fab model, we refer the reader to [11].
For the purpose of demonstration, we replace the observed system by a discrete event simulation. We take an influx profile that we get out of the averaged DES models and then we try to predict the non-equilibrium behavior using the PDE model. In the simulations, we compute the transport coefficients depending on a macroscopic state variable (Z), which is the total WIP since we use FIFO (first in first out) as the service rule in the system. One can also use PULL or PUSH policies as service rules; then it would be wise to use upstream WIP or downstream WIP in the simulations.
We execute the DES model by generating a total of 900 lots with a varying influx. In Figure 1, the averaged influx computed out of 200 discrete event simulations is seen. We compute the data for the observations from a steady state situation and then predict the system in a non-steady state regime as can be seen by the transients in Figures 4 and 5.
Velocity coefficients for the model of semiconductor fab is shown in Figure 2 both in contour and log scale plots as functions of WIP for each stage. Moreover, diffusion coefficients is shown in Figure 3 as functions of WIP for each stage in contour and log scale plots. After estimating the velocity and the diffusion coefficients, which are extracted from the 'observed data' as outlined in Section 3, we numerically solve the PDE given in (1) using an upwinding scheme where the flux is given as F = Cρ − D∂ x ρ. With the PDE model, we predict the transient behavior of the DES model according to the influx given in Figure 1 that we get out of the average of 100 DES model runs. The total time to execute this DES model takes almost one day, whereas the time to simulate the corresponding PDE model takes around one hour. The execution in this example is carried out by MATLAB 5.01 and obtained by a PC Pentium 4, 2.00 6 Hz, 1 GB of RAM.
The comparisons of all DES and PDE fluxes as functions of time at each machine stage can be seen in 3D surface plots in Figure 4. Similarly, the WIP levels of DES and PDE models are compared as functions of time at each machine stage in 3D surface plots in Figure 5. In the simulations, each time step corresponds to 120 mesh units depending on what time unit is chosen; i.e. minutes, hours, or days, etc.
As it can be seen in Figures 4 and 5, we substantiate the agreement of PDE and DES models by comparing the WIP levels and flux values of the two models at each stage. In addition, at two randomly chosen sample stages (8 and 23), we compare three averaged discrete event experiments to the PDE model in 2D plots in terms of WIP levels as functions of time. It is seen in Figures 6 and 7 that, in each plot, all the curves fall in the same imaginary envelope. Therefore, using the PDE model we are able to predict the transient behavior of the given system with the macroscopic model under conditions different than those used to collect the data.  Conclusion. In this paper, we have presented a computationally efficient way to simulate stochastic production systems. Our goal was to solve the diffusive PDE model (1) after extracting the transport coefficients (C and D) from the observed data which we obtained out of a discrete event experiments for the purpose of this paper. With the PDE model, using a MATLAB program, we predicted the transient behavior of the DES model according to the averaged influx obtained out of 100 discrete event experiments. Finally, we compared the PDE and DES models  in terms of WIP levels and fluxes at each machine stage as functions of time and substantiated the agreement of two models. PDE models are preferable to simulate large and complex production systems, because they are accurate, approximately twenty times faster than DES models in execution of simulations, and they can be optimized via constraint optimization, which is a very important advantage for a mathematical model as compared to other models.