A HYPERBOLIC RELAXATION MODEL FOR PRODUCT FLOW IN COMPLEX PRODUCTION NETWORKS

This paper presents 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, in a relaxation limit, exhibit the correct diffusive properties 
given by the variance of the observed data.

1. Introduction. We present a continuum -traffic flow like -model for the evolution of a large number of parts through a complex production system. The basic idea of traffic flow like models for production systems is to model the evolution of parts as moving on a virtual 'road', introducing an artificial spatial variable, namely the stage of the production process. So, the raw part enters the system at stage x = 0 and the finished part leaves the system at stage x = 1. In the work presented in this paper we compute the velocity of the part from the statistics of a large number of observed data. This leads to the formulation of a Lagrangian particle model, and the consequent derivation of macroscopic conservation laws for the part density using large time averages. The present paper represents an extension of the work in [1] and [5]. In this previous work, the large time averaged conservation laws were derived from a standard Chapman -Enskog expansion of an underlying kinetic equation, leading to a single parabolic conservation law, i.e. a convection diffusion equation. As will be explained in the following, the parabolicity of the resulting equation has to be viewed as an artifact of the asymptotics and actually gives rise to modeling problems, as far as the formulation of boundary conditions is concerned. In the present paper, the parabolic conservation law is replaced by a hyperbolic relaxation system which, on one hand, exhibits the same unidirectional flow of information as the underlying particle model, and, on the other hand, reproduces the parabolic equations derived and discussed in [1] and [5] in the limit of very large systems.
This paper is organized as follows: In Section 2, we explain in detail how observed data are used to compute state dependent probability distributions for the velocities. Section 3 represents a review of the derivation of the particle model and the resulting kinetic equation for the particle distribution function in [1]. Section 4 is devoted to the derivation of the hyperbolic relaxation model and its relation to the parabolic model discussed in [1] and [5]. Some numerical results on a suitably simplified model of a semiconductor fab are presented in Section 5.
2. Parameter Extraction. The goal of this paper is to develop a macroscopic model for the flow of individual parts through the production process. Other than in the standard discrete event simulation (DES) [2] models, we will not try to model the whole process by a synthesis of -necessarily simplified -models for individual nodes, but instead present a methodology for extracting the transport coefficients from actually observed data. We assume, that we have observed the evolution of a large number of parts through the system over a period of time and have recorded the times each part has entered each stage of the process. So, the data available to the model are in the form of a matrix {a k (n), k = 1 : K + 1, n = 1 : N }, where a k (n), k = 1 : K denotes the time part number n has entered stage number k of the process, and a K+1 (n) denotes the time the (finished) part number n exits the whole system. We will employ a traffic flow type model. That is, we will model the progress of each part through the system as the evolution of a particle on a virtual 'road', located in the interval [0, 1]. So parts enter at stage x = 0 and leave at stage The goal here is to develop a non -equilibrium model, i.e. a model which is capable of responding to transient inputs. We assume that this response is essentially given by the response to some integral quantities S k . So, at the same time we observe the -possibly vector valued -quantity s k (n) ∈ R d which denotes the measurement of S k at the time a k (n). We employ a mean field assumption, that is, we assume that the evolution of an individual part does not influence the quantity S k significantly, and therefore the random variables S k are in zeroth order independent of the observed velocities V k . This yields a joint probability distribution for v k and S k of the form In the following, we will compute the macroscopic states S k from a homogenized density for the individual parts. The equation for this homogenized density will employ the conditional probability dP[V k = v | S k = s] of the velocity V k , given that the state S k . As given by (1), this conditional probability is not well defined, since the joint probability density dP[V k = v, S k = s] is discrete in the variable s, i.e. the probability to observe a given state S k is either zero or infinity. We therefore modify the definition (1) by distributing the observed variables S k over cells in R d .
Let R d be the the disjoint union of cells C j , j ∈ N. We replace (1) by where χ j denotes the indicator function of the cell C j ⊂ R d , and |C j | = χ j (s) ds denotes the volume of the cell C j . The definition (2) allows us to define the conditional probability density of the velocity V k , given the observation S k in the usual way as dP dv is a probability density in the variable v, which depends on the macroscopic state density S(x, t). How well the probability distribution U S in (4) is able to predict the evolution of parts depends, of course, on the actual degree of correlation between the velocities v k and the observed states s k . The appropriate choice of the observed states s k is therefore essential. In practice, we will choose the s k as either the total load (the total number of parts in the system) or as the number of parts at a stage ≥ k or ≤ k, or a combination of the above (see [5] for details). We note, that U s (x, v, t) represents in some sense a generalization of the concept of clearing functions [4] to clearing distributions.
3. The Particle Model. This section is essentially a summary of the work in [1], leading to a kinetic equation for the evolution of the density f (x, v, t) of parts at stage x having velocity v. Denoting the position of the part at time t by x = ξ(t), the evolution of the part is given by the Newton equation dξ dt = v where the velocity v is updated periodically from the distribution U S defined in Section 2. This implies that the position ξ after the infinitesimal time interval dt is given by We define the frequency with which these updates are performed by ω and obtain the following rule for computing the velocity at time t + dt from the velocity at time t. At each infinitesimal time interval dt we 'flip a coin', that is we compute a random variable r, where the probability that r = 1 equals ωdt and the probability that r = 0 equals 1 − ωdt. If r = 0, we maintain the current velocity and, if r = 1, we update the velocity v from the distribution U defined in Section 2. This gives the rule (5) and (6) define a stochastic process for the evolution of the particle position ξ(t) and its velocity v(t). If we define the probability density that the part at time t is at position x with velocity v by f (x, v, t) dxdv, we obtain the evolution equation In the limit dt → 0 this gives the evolution equation (see [1] for details) Note that we have made the scattering frequency ω, defined in (6), dependent on the current velocity v as well as on the position x. The reason for this dependence is that the scattering velocity ω has to be defined in terms of a mean free path, i.e. the distance λ a part travels before it undergoes a change in velocity. If we update, on average, each time a part enters the next cell, we have λ = ε. Since the average time between a sudden change in velocity is, according to (6), given by 1 ω , a part with velocity v travels, on average, a distance λ = v ω before undergoing a scattering event. This gives ω = v ε , and the collision operator Q in (7) becomes We assume a complex system with a large number K of stages, and so ε = 1 K << 1 holds. The goal is therefore to derive an evolution equation for the part density ρ(x, t) = f (x, v, t) dv, using some form of asymptotics for ε → 0. Note that the propagation of information in the kinetic equation (7) is unidirectional, that is parts can only move forward and the support of the density f (x, v, t) is always confined to the half space v > 0, since velocity can only change randomly according to U S whose support is also confined to v > 0. It therefore suffices to prescribe the influx density f (x = 0, v, t) for positive velocities v to obtain a well posed boundary value problem for the kinetic equation (7).

4.
Asymptotics. The goal of this section is to derive a simple conservation law for the density of parts, given by ρ(x, t) = f (x, v, t) dv. In [1], this is achieved by a Chapman -Enskog expansion of the kinetic equation (7). (We refer the reader to [3] for background on Chapman -Enskog asymptotics.) This leads to a parabolic conservation law for the density ρ. The parabolicity of the macroscopic model represents, however a major shortfall. In a parabolic model information is transmitted at infinite velocity in both spatial directions. This is not the case for the underlying particle model in Section 3. Amongst other things, this implies that the parabolic model will require a boundary condition at the outflux boundary x = 1, and there is no physical way to prescribe the outflux. In [1] and [5] this problem has been dealt with by using artificial numerical boundary conditions. The deeper reason for this conundrum is that Chapman -Enskog asymptotics translates the diffusion in velocities, produced by the collision operator Q in (7) into diffusion in the spatial variable x, and is , in a sense, too coarse an instrument to reflect that this diffusion only takes place among positive velocities. We will therefore employ a modified Chapman -Enskog procedure, leading to a hyperbolic system with positive characteristic velocities, which in a relaxation limit reduces to the parabolic system in [1].
The basic idea is to split the dynamics of equation (7) into a slow part, constrained to the kernel of the collision operator Q, and a fast part constrained to the complement of the kernel. We first characterize the kernel of the collision operator Q in (8). The kernel of Q is given by functions of the form where c is some arbitrary function of space and time. For the following it will be convenient to define the normalized kernel distribution as We split the function space B for the kinetic density function f into the linear hull of the kernel and its complement: We define the projection onto the kernel space B 0 by The definition (9) of the projection operator P implies immediately that P Q = QP = 0 holds. We split the evolution equation (7) into the evolution in the kernel manifold and the evolution in its complement B 1 by writing , and by applying the projections P and I − P to equation (7), giving Equation (10)(a) gives the conservation law for the density ρ, whereas equation (10)(b) has to be solved asymptotically for ψ, giving the closure for the conservation law (10)(a). In detail we have A closed conservation law for the density ρ is obtained by computing the term vψ dv in (11)(b) asymptotically from the equation (11)(d). In the usual Chapman -Enskog procedure, this is done by using the pseudo -inverse Q + of the operator Q, i.e. the inverse of Q restricted to B 1 . A direct calculation gives Q + [g] = (I − P )[ g v ], and (11)(d) can be written as A. UNVER, C. RINGHOFER AND D.ARMBRUSTER which yields, in first order Inserting this into the closure in (11)(b) gives (13) and the asymptotic conservation law The convection diffusion equation (14) is the the parabolic conservation law derived in [1]. Computing the diffusion coefficient, using the definition (9) of the projection P , gives D = w 1 (w 1 w −1 − 1), where from now on we use the notation w j = v j W dv, j ∈ Z. The diffusion coefficient D can be written in terms of a variance as and is therefore strictly positive. The conservation law (14) has been used in [1] and [5] to model quite complex, re -entrant production systems. It does, however, exhibit the fundamental flaw outlined in the beginning of this section. In the present work, we will take a different approach, replacing the diffusion approximation (14) by a hyperbolic relaxation model. Instead of approximating the solution of (11)(d) using the pseudo inverse of the operator Q, we approximate the component ψ of the density f , belonging to B 1 by a moment closure. So, we approximate ψ(x, v, t) by the ansatz ψ(x, v, t) = β(x, t)R(x, v, t) for a given density R and replace (11)(d) by its first moment, giving We assume R to be normalized, such that vR dv = 1 holds. This gives the equation

and the hyperbolic system
with the system matrix A and the coefficient q given by The hyperbolic system (16) will depend of course on the choice of the shape function R. In order to match the hyperbolic system as closely as possible to the original Chapman -Enskog solution (12), we choose the shape function R to be in highest order equal to the kernel element (12) given by the Chapman -Enskog expansion. Therefore we set where the moments of the equilibrium density W are defined as w j (x, t) = v j W (x, v, t) dv and D denotes the diffusion coefficient given in (15).
Proof. Computing the projections, we obtain and the normalization coefficient γ is therefore given by with D the diffusion coefficient, defined in (15). This gives (19).
As it turns out, this choice of the shape function R yields a system matrix A in the hyperbolic relaxation system (16) which retains the property of the original particle system, namely that all characteristic velocities are non -negative, i.e. that both eigenvalues of the matrix A in (16) are nonnegative. To this end, it is convenient to introduce the variance σ 2 of the distribution W . We define We have Lemma 4.2. The system matrix A and the coefficient q in (17) is given by Both eigenvalues of A are real and nonnegative.
Proof. Computing the (2, 1) component of the matrix A in (17), we obtain Computing the (2, 2) term of the matrix A in (17) we obtain, using the definition To compute the coefficient q, we have, according to (19) and the definition of the collision operator Q,

5.
A Computational Example. We carry out a numerical experiment in order to compare a complex system to the hyperbolic relaxation 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. We refer the reader to [5] for a detailed explanation of the simplified semiconductor fab model.
For the purpose of demonstration, we replace the observed system by a discrete event simulation (DES). 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 hyperbolic relaxation model.
We execute the DES model by generating a total of 900 lots with a varying influx. The averaged influx, computed out of 100 discrete event simulations, can be seen in Figure 1. Afterwards, we compute the data for the observations from various steady state situations and then try to predict the system in a non -steady state regime as can be seen by the transients in Figures 2 and 3.
As explained in section 4, after computing the transport coefficients in the matrix A and the coefficient q, given simply in (21), we solve the hyperbolic relaxation system (16) by using upwinding scheme.
As we see in Figures 2 and 3, although the WIP levels and flux values of hyperbolic relaxation and DES models at each stage are not exactly equal, the results are close enough to each other to substantiate the agreement of two models. Therefore, 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.