An Upwind-Euler Scheme for an ODE-PDE Model of Supply Chains

An Upwind-Euler scheme for an ODE-PDE model of supply chains (1) Alfredo Cutolo (1) , Benedetto Piccoli (2) , Luigi Rarità (1) Department of Information Engineering and Applied Mathematics, University of Salerno, Fisciano (SA), Italy (2) Istituto per le Applicazioni del Calcolo “Mauro Picone”, Consiglio Nazionale delle Ricerche, Roma, Italy August 30, 2010 Abstract Numerical techniques for the simulation of an ODE-PDE model for supply chains are presented. First we describe a scheme, based on Upwind and explicit Euler methods, provide corrections to maintain positivity of solutions, prove convergence and provide convergence rate. The latter is achieved via comparison with Wave Front Tracking solutions and the use of generalized tangent vectors. Different choice of time and space meshes give similar results, both for CPU times and numerical errors. Fast algorithms, based on an accurate choice of time and space meshes and data structures, are then proposed, achieving high computational gains. 1 Introduction In last years, scientific communities showed a great interest in modelling dynamics of industrial productions, managed by supply chains. This study is fundamental in order to reduce some unwished phenomena (bottlenecks, dead times, and so on). Several mathematical approaches have been proposed. For example, some models are discrete and based on considerations of individual parts (for a review, see [10]). Other models are continuous ([1, 2, 3, 24]), and some based on partial differential equations. The first paper in this last direction was [1], where the authors, via a limit procedure on the number of parts and suppliers, have obtained a conservation law formulation. The flux involves parts density and the maximal productive capacity. The main problem is that solutions for discontinuous productive capacity functions may exhibit delta waves. For this reason, alternative models were proposed. For instance, a mixed continuum-discrete model is given in [11]: The supply chain is described by continuous arcs and discrete nodes. Also extensions on networks have been made ([12, 13, 14]). In this paper, we focus the attention on the continuous model for supply chains and networks proposed by Göttlich, Herty and Klar in [17]. Our main aim is to analyze the convergence of the Upwind-Euler scheme and improve its performance by accurate choice of discretization parameters. Examples of schemes for conservation laws on networks can 1 be found in [6, 17, 20]. The necessity of having fast numerics is justified by optimization problems which arise naturally in applications, see [18, 22]. The model of [17] is obtained considering suppliers on which the processing rate is constant (thus avoiding the problem of delta waves) and having queues in front of each supplier. The evolution of the queue buffer occupancy is given by the difference of fluxes from the preceding and following suppliers. The outcome is a coupled ODE-PDE model. Two numerical schemes must be considered: One for the ODE and one for the PDE. We choose the Upwind method for parts densities, described by conservation laws (PDEs), and the explicit Euler scheme for queue evolutions, modelled by ODEs. For details, see also [16] and [23]. Notice that the ODE-PDE model guarantees positivity of queue buffer occupancies (and densities), while this is not granted for the Upwind-Euler scheme. Thus first we consider fluxes corrections to avoid this drawback. The choice of time and space meshes can be uniform over suppliers only in case of rational ratios among lengths. Thus we consider different discretizations of the Upwind-Euler to deal with the general case and also to reduce computational complexity. Convergence of the scheme is proved using a comparison with Wave Front Tracking approximate solutions. More precisely, we consider nesting grids in which meshes are divided by two in successive approximations. Then generalized tangent vectors, as in [21], are defined and the evolution of their norms permits to bound L1 distances of solutions. Finally, a linear convergence rate is achieved in terms of the space (or time) mesh. A similar technique was used in [8] to prove convergence for an ODE-PDE (not completely coupled) model to track the trajectory of a car on a road network. From several simulations of supply chains, it is evident that heavy numerical errors occur if corrections for negative queues are not adopted. Moreover, different time-space discretizations for the Upwind-Euler method perform similarly, both in terms of CPU times and numerical errors. Focusing on the last point, computational gains are obtained by the use of logic pointers to update the computed density values. This allows a strong reduction of computational complexity as shown by simulation of relatively large supply chains (10 and 100 arcs), for which CPU times reduce of up to 90%. The outline of the paper is the following. In Section 2, we present the ODE-PDE model. Then Section 3 describes the Upwind-Euler scheme, while its convergence is proven in Section 4, shifting to the Appendix the description of generalized tangent vectors. Section 5 reports tests showing the importance of positivity of solutions. Reduction of computational complexity is addressed in Section 6 and the relative numerical tests are described in Section 7. 2 A model for supply chains In this section, we present an ODE-PDE model for supply chains first proposed in [17], and based on the work [1]. Beside the conservation laws formulation proposed in [1], such model includes time - dependent queues describing the transition of parts among suppliers. A supply chain is a directed graph consisting of arcs J = {1, ..., N } and vertices 2 V = {1, ..., N − 1}. Each arc j ∈ J , parameterized by an interval [aj , bj ], models a supplier. Each vertex is connected to one incoming arc and one outgoing arc and we assume that arcs are consecutively labelled, i.e. arc j is connected to arc j +1 and bj = aj+1 (see Figure 1). For the first and the last arc, we either set a1 = −∞ and bN = +∞, respectively, or provide boundary data for the inflow and outflow. Each supplier j has a queue in front of itself, i.e. at x = aj , and is characterized by a qj j 1 b1 b j 1 a2 aj bj a j 1 Figure 1: an example of supply chain. maximum processing capacity µj > 0, length Lj > 0 and a processing time Tj > 0. The rate Lj /Tj represents the processing velocity. Indicating by ρj (t, x) the density of parts in the supplier j at point x and time t, the evolution of parts is described by a conservation law: ∂t ρj (t, x) + ∂x fj (ρj (t, x)) = 0, ∀ x ∈ [aj , bj ] , t > 0, (1) where the flux function fj (ρj (t, x)) is given by: fj : [0, +∞[ → [0, µj ] , ¾ ½ Lj fj (ρj (t, x)) = min µj , ρj (t, x) . Tj The interpretation of equation (1) is the following: Parts are processed with velocity Lj /Tj but with a maximal flux µj . Each queue buffer occupancy is a time dependent function t → qj (t). If the capacity of the supplier j − 1 and the demand of the supplier j are not equal, the queue increases or decreases its buffer. More precisely, we have: d qj (t) = fj−1 (ρj−1 (bj−1 , t)) − fj (ρj (aj , t)) , dt j = 2, ..., N, (2) qj (0) = qj,0 ≥ 0, where fj−1 (ρj−1 (bj−1 , t)) is defined by the evolution on supplier j − 1, while the flux on the outgoing arc j is defined as ½ min{fj−1 (ρj−1 (bj−1 , t)) , µj } , qj (t) = 0, fj (ρj (aj , t)) = (3) µj , qj (t) > 0. Notice that the flux fj (ρj (aj , t)) depends on the capacity of the queue: If the queue buffer is empty, the inflow to supplier j is equal to the outflow from supplier j − 1, otherwise the inflow is maximal. In other words, when there are parts in the queue the processing 3 occurs at the maximal possible rate, namely µj , so to empty the queue as fast as possible. Finally, the complete system of equations is: ½ ¾ Lj ∂t ρj (t, x) + ∂x min µj , ρj (t, x) = 0, ∀ x ∈ [aj , bj ] , t > 0, j ∈ J , (4) Tj ρj (0, x) = ρj,0 (x) ≥ 0, ∀ x ∈ [aj , bj ] , d qj (t) = fj−1 (ρj−1 (bj−1 , t)) − fj (ρj (aj , t)) , j = 2, ..., N, dt qj (0) = qj,0 ≥ 0, ½ min{fj−1 (ρj−1 (bj−1 , t)) , µj } , qj (t) = 0, fj (ρj (aj , t)) = µj , qj (t) > 0. (5) (6) (7) (8) Lemma 1 Consider a supply chain evolution ρj (t, x), qj (t), i.e. a solution to (4)-(8). Then for every j ∈ J , t ≥ 0 and x, it holds ρj (t, x) ≥ 0, qj (t) ≥ 0. Proof. Since the inflows (8) are positive and the fluxes fj vanish at 0, the density parts ρj are always positive by comparison principle of conservation laws. Moreover (6) and (8) guarantee that the derivative of queue buffer occupancy is always positive when the queue is empty, thus the conclusion follows. 2.1 Supply networks. The ODE-PDE model (4)-(8) can be extended to the case of networks, see [18]. The main idea is to introduce traffic distribution coefficients at nodes, which describe how the outgoing flux from a given node distribute over suppliers which are downstream. More precisely, a network is a directed graph formed by a set of arcs J and a set of vertices V. Each arc j ∈ J , parameterized by an interval [aj , bj ], models a supplier (possibly having infinite endpoints) and there is a queue in front of it if aj > −∞. Each vertex v is connected to some incoming arcs and some outgoing arcs. On the other side each arc is outgoing for at most one vertex and is incoming for at most one vertex. More precisely j is outgoing for some vertex if aj > −∞ and is incoming for some vertex if bj < +∞. We indicate by Inc(v) ⊂ J the set of incoming arcs into v and by Out(v) ⊂ J the set of outgoing arcs from P v. For each v there exist distributions coefficients (av,j )j∈Out(v) such that av,j ∈]0, 1[ and j∈Out(v) aj,v = 1. The coefficient av,j represents the percentage of flux outgoing from v which is directed to the supplier j. The model (4)-(8) can be modified as follows to deal with networks. Notice that (4), (5) and (7) have the same expression. Now fix a supplier j, with aj > −∞, and let v be the vertex such that j ∈ Out(v). Equation (6) is replaced by: d qj (t) = av,j fv (t) − fj (ρj (aj , t)) dt where fv (t) = X fk (ρk (bk , t)) k∈Inc(v) 4 (9) (10) and (8) is replaced by: ½ fj (ρj (aj , t)) = min{av,j fv (t), µj } , µj , qj (t) = 0, qj (t) > 0. (11) Also in this case we have positivity of solutions, more precisely: Lemma 2 Consider a supply network evolution ρj (t, x), qj (t), i.e. a solution to (4), (5), (7), (9), (10), (11). Then for every j ∈ J , t ≥ 0 and x, it holds ρj (t, x) ≥ 0, qj (t) ≥ 0. Proof. Since aj,v > 0, the inflows (11) are positive and, as before, the fluxes fj vanish at 0, thus the density parts ρj are always positive by comparison principle of conservation laws. Equations (9), (10) and (11) ensure that the derivative of queue buffer occupancy is always positive when the queue is empty, thus the conclusion follows. 3 The Upwind-Euler scheme In this Section we introduce the Upwind-Euler method with various possible discretizations. For simplicity we focus on supply chains, being the case of networks entirely similar. Numerical results for parts dynamics on a supply chain are obtained finding, for each arc j, suitable approximation for the density ρj (t, x), and the queue buffer occupancy qj (t). In particular, we need a PDE numerical method and an ODE one: We choose the Upwind scheme for the PDE and the explicit Euler scheme for the ODE. For each arc j ∈ J , define a numerical grid in [0, Lj ] × [0, T ] using the following notations: L • ∆xj = Njj is the space grid mesh, where Nj is the number of segments into which we divide the j-th supplier; • ∆tj = ηTj is the time grid mesh, where ηj denotes the number of segments into which [0, T ] is divided; • (xi , tn )j = (i∆xj , n∆tj ), i = 0, ..., Nj , n = 0, ..., ηj are the grid points; • j ρni is the value taken by the approximated density at the point (xi , tn )j ; • qjn is the value taken by the approximate queue buffer occupancy at time tn . The Upwind method reads: j n+1 ρi = j n ρi − ¢ ∆t Lj ¡ j n ρi − j ρni−1 , j ∈ J , i = 0, ..., Nj , n = 0, ..., ηj , ∆x Tj with CFL condition given by ∆t ≤ ∆x maxj Lj Tj , while the explicit Euler method is given by: ¡ n ¢ n qjn+1 = qjn + ∆t fj−1,out − fj,inc , j ∈ J \ {1} , n = 0, ..., ηj , 5 (12) (13) (14) where n fj,inc = n fj−1,out = fj−1 (j−1 ρnNj−1 ), n o ( min fj−1 (j−1 ρnNj−1 ), µj , qjn (t) = 0, qjn (t) > 0. µj , (15) Boundary data are treated using ghost cells and the expression of inflows given by (15). Assuming that the numbers Lj have rational ratios, it is possible to choose a space grid mesh ∆x and a common time grid mesh ∆t. We will discuss below the general case. 3.1 Correction of numerical fluxes in case of negative queues The ODE numerical scheme does not necessarily maintain the positivity properties of the true solutions given by Lemma 1. We thus modify the Euler scheme so as to accomplish positivity of queue buffer occupancies. £ £ Consider a supplier j and a time interval tn , tn+1 so that qjn+1 < 0 < qjn . Then we define qj (t) for every time t by linear interpolation, see Figure 2, namely qj (t) = qjn+1 − qjn ∆t t+ qjn tn+1 − qjn+1 tn ∆t , £ £ t ∈ tn , tn+1 . (16) Then q(·) vanishes at some time t > tn , which is computed by (16): qj t 't ' qj n t n 1 tn t t q j n 1 't Figure 2: Negative queue buffer occupancy at tn+1 . t = tn + ∆t0 , ∆t0 = qjn n+1 ∆t = qjn − qj qjn . n µj − fj−1,out (17) £ ¤ Correcting to zero qj (t) , t ∈ t, tn+1 , the following numerical correction for the entering n flux fj,inc is needed: n fj,inc = ¡ ¢ n ¤ 1 £ 0 ∆t µj + ∆t − ∆t0 fj−1,out . ∆t 6 (18) The correction (18) on the boundary incoming data for the supplier j influences the approximation of ρj (x, t), with consequent effects on the dynamics for following suppliers and queues. We have the following (cfr. Lemma 1): Lemma 3 Consider a numerical solution j ρni , qjn , defined using the flux correction (18), such that j ρ0i ≥ 0 and qj0 ≥ for all j ∈ J , i = 1, . . . , Nj . Then j ρni ≥ 0, qjn ≥ 0, for all n ≥ 0, j ∈ J , i = 1, . . . , Nj . Proof. By formulas (12), (13) and (15), since j ρ0i ≥ 0, for all j ∈ J , i = 1, . . . , Nj , we clearly have that j ρni ≥ 0, qjn ≥ 0, for all n ≥ 0, j ∈ J , i = 1, . . . , Nj . Now if (14) gives rise to a negative value for some j ∈ J and n ≥ 1, then using the flux correction (18), we can write: µ ¶ ¡ n ¢ ¡ ¢ n ¤ 1 £ 0 n+1 n n n n 0 q = qj +∆t fj−1,out − fj,inc = qj +∆t fj−1,out − ∆t µj + ∆t − ∆t fj−1,out . ∆t Then by (17), we have: £ ¡ ¢ n ¤ n n q n+1 = qjn + ∆tfj−1,out − ∆t0 µj + ∆t − ∆t0 fj−1,out = qjn + ∆t0 (fj−1,out − µj ) = 0. Remark 4 An alternative method to avoid negativity of queues is the use of adaptive time meshes, where ∆t is replaced by ∆t0 computed in (17). Alternatively, one may modify the equation for q with a relaxation term when the queue buffer is decreasing to zero. In other words q is exponentially decaying to zero when the queue is emptying. This gives rise to a nonlinear stiff equation and restrictions on the time steps must be carefully addressed, see [19]. 3.2 Different space and time grid meshes We now consider the possibility of choosing different space and/or time grid meshes. This is necessary in the general case in which Lj have not rational ratios, but can also be useful for computational complexity reduction, as we will see in Section 6. 3.2.1 Different space meshes for different suppliers For each supplier j ∈ J , the numerical grid in [0, Lj ] × [0, T ] is defined choosing a fixed time grid mesh ∆t, then different space grid meshes are necessary. Set ∆xj = vj ∆t, where L vj := Tjj indicates the processing velocity. Then, grid points are (xi , tn )j = (i∆xj , n∆t), i = 0, ..., Nj , n = 0, ..., ηj . The Upwind scheme for the parts density of the arc j now reads: j n+1 ρi = j n ρi − ¢ ∆t ¡j n vj ρi − j ρni−1 , j ∈ J , i = 0, ..., Nj , n = 0, ..., ηj . ∆xj 7 (19) The CFL condition (see [23]) is automatically satisfied since: ½ ¾ ∆xj ∆t = min :j∈J . vj (20) For queues we refer again to equation (14). In case of negative values of queues, fluxes corrections have to be considered, which are the same as in Section 3.1. 3.2.2 Different time meshes for different suppliers We now consider the case of different temporal meshes (and same space meshes). Fix two consecutive arcs j − 1 and j, then two different numerical grids are defined, whose points are, respectively: (xk , tnj−1 )j−1 = (k∆x, nj−1 ∆tj−1 ) , k = 0, .., Nj−1 , nj−1 = 0, ..., ηj−1 , ¡j ¢ xk , tnj = (k∆x, nj ∆tj ) , k = 0, .., Nj , nj = 0, ..., ηj . For the queue buffer occupancy the explicit Euler is given by: ³ ´ n +1 n nj nj qj j = qj j + ∆tj fj−1,out − fj,inc . n (21) n j j where fj,inc is computed as in (15), while fj−1,out must be suitably defined. If ∆tj−1 < ∆tj , we define m (nj ) and M (nj ) as: m (nj ) = sup {m : m∆tj−1 ≤ nj ∆tj } ; M (nj ) = inf {M : M ∆tj−1 ≥ (nj + 1) ∆tj } , and set n j fj−1,out = M (nj )−m(nj )−1 X m(n )+l ∆tj−1 fj−1 (j−1 ρNj−1j ) l=1 m(n ) + [(m (nj ) + 1) ∆tj−1 − nj ∆tj ] fj−1 (j−1 ρNj−1j )+ M (n )−1 + [(nj + 1) ∆tj − (M (nj ) − 1) ∆tj−1 ] fj−1 (j−1 ρNj−1j ). Notice that, in the special case ∆tj = γ∆tj−1 , γ ∈ N r {1}, we simply have: n j fj−1,out = M (nj )−m(nj )−1 X m(n )+l ∆tj−1 fj−1 (j−1 ρNj−1j )= l=1 γ X γn +l j ∆tj−1 fj−1 (j−1 ρNj−1 ). l=1 If, on the contrary, ∆tj−1 > ∆tj , we set  n nj ∆tj ∆tj−1 j fj−1,out = fj−1 where b·c indicates the floor function. CFL conditions now reads: ∆x ∆tj ≤ , vj  , j = 1, . . . , N. In case of negative values of queues, fluxes corrections have to be considered. The modifications w.r.t. Section 3.1 can be easily computed and are omitted for sake of space. 8 4 Convergence The aim of this Section is to study the convergence of the Upwind-Euler numerical scheme. The main idea is to compare numerical solutions with those produced by a Wave Front Tracking algorithm and control the norm of generalized tangent vectors which measure the distance, as in [21]. We consider a Cauchy problem of the type (4)-(8), with initial conditions ρj,0 in BV, the space of bounded variation functions. Fix an initial space meshes ∆x0 and define a sequence of approximate solutions ν,j ρni , generated sampling the initial datum ρj,0 on grids with space meshes ∆xj,ν = 2−ν ∆xj,0 and using the time meshes: ∆tj,ν = ∆xj,ν ∆xj,0 = 2−ν , vj vj (22) where vj is the velocity of the j-th suppliers, thus granting the CFL conditions. More precisely: ¡ ¢ ν,j 0 ρi = ρj,0 (aj + i 2−ν ∆xj,0 ) + where (· +) indicates the limit from the right, which exists because of the assumption of BV initial data. We can define a projection of the approximate solution over the space of piecewise constant functions by setting: πP C ¡ν,j n ρ ¢ Lj /2−ν ∆xj,0 −1 = X ν,j n ρi χ[aj +i 2−ν ∆xj,0 ,aj +(i+1)2−ν ∆xj,0 [ i=0 where χ[a,b] is the indicator function of the set [a, b]. We define the corresponding queue buffer occupancy approximations ν qjn as specified in the previous Section and consider the projection over piecewise linear function: πP L (ν qj ) (t) = ν qjn + (t − n∆tj,ν )(ν qjn+1 − ν qjn ) for t ∈ [n∆tj,ν , (n + 1)∆tj,ν [. F T starting We will also consider the Wave Front Tracking (briefly WFT) solution j ρW ν from the initial datum ¡ ¢ πP C ν,j ρ0 . Roughly speaking, a WFT solution is obtained in the following way: • Solve the Riemann problems corresponding to discontinuities of πP C ¡ν,j ¢ ρ0i ; • Use the piecewise constant solution obtained piecing together the solutions to Riemann problems up to the first time of interaction of two shocks; • Then solve a new Riemann problem created by interaction of waves and prolong the solution up to next interaction time, and so on. To ensure the existence of WFT solutions and their convergence, it is enough to control the number of waves, interactions and the BV norm. This is easily done in scalar case since both the number of waves and the BV norm are decreasing in time. We refer the 9 reader to [4] for details. For queues we use the exact solutions to (6) and indicate them by ν qjW F T . BV estimates for complete ODE-PDE model (4)-(8) are proved in [21]. For a general WFT scheme we should replace, in solutions to Riemann problems, possible rarefactions by a set of small non entropic shocks of size 2−ν . However, in our case rarefactions do not show up because the flux is piecewise linear. Let us describe in detail solutions to Riemann problems. Fix a supplier j and consider a Riemann problem with initial data (ρ− , ρ+ ). Then we have to distinguish some cases: • ρ− < µj . In this case the solution is given by a shock travelling with velocity f (ρ )−f (ρ ) λ = j ρ++ −ρj− − (which equals vj in case ρ+ ≤ µj ). • ρ− = µj . Then the solution is a shock travelling with velocity vj if ρ+ < µj and with zero velocity if ρ+ > µj . • ρ− > µj . If ρ+ ≥ µj then the solution is a shock with zero velocity. Otherwise, the solution is formed by a first shock (ρ− , µj ) travelling with zero velocity and a second shock (µj , ρ+ ) travelling with velocity vj . Notice that as soon as a boundary datum will achieve a value below µj , then in finite time all values above µj will disappear from the j-th supplier, see also [21]. Therefore, for simplicity, we will assume ρj,0 (x) ≤ µj . (23) Then the same inequality will be satisfied for all times both for the numerical approximation and wave front tracking ones, i.e.: ¡ ¢ πP C ν,j ρn (x) ≤ µj for all n ≥ 0, x ∈ [aj , bj ], (24) j WFT ρν (t, x) ≤ µj for all t ≥ 0, x ∈ [aj , bj ]. (25) In this case solutions to Riemann problems are particularly simple, indeed the conservation law is linear, thus given some Riemann data (ρ− , ρ+ ) on the j-th supplier, the solution is always given by a shock travelling with velocity vj . Therefore, the WFT solution will differ from the projected numerical ones only for effects due to queues dynamics. To prove convergence numeri¡ ¢ of the numerical scheme, we thus compare the projected F T ,ν q W F T ). cal solution (πP C ν,j ρn , πP L (ν qj )) with the wave front tracking solution (j ρW ν j More precisely we consider their difference as expressed by a generalized tangent vector, i.e. write: X ¡ ¢ FT πP C ν,j ρn (t) = j ρW (t) + ξα (t)∆ρα (t), (26) ν α∈Aj (t) πP L (ν qj )(t) = ν qjW F T (t) + ηj (t), (27) F T at time t and ∆ρ the strength of where Aj (t) is the set of discontinuities of j ρW α ν the given¡ discontinuity (i.e. the value of the jump). In other words the part den¢ ν,j n j W F T sity πP C (t) shifting the shocks of quantities ξα , while ρ (t) is obtained from ρν πP L (ν qj )(t) is obtained from ν qjW F T shifting its value of ηj (t). We use the symbol (ξ, η)ν 10 to indicate the collection of all tangent vectors over all suppliers. This approach was used in [21] to study uniqueness and continuous dependence of solutions. We report the general construction in the Appendix for reader’s convenience. The norms of tangent vectors to WFT solutions are proved to be decreasing in time as stated in Lemma 8 of the Appendix (for a proof see [21]). Our idea is to study the evolution of tangent vectors defined in (26)-(27) to estimate the increase in time of their norms defined as: X X X |ηj (t)|. |ξα ∆ρα | + k(ξ, η)ν (t)k = j j α∈Aj (t) Notice that: k(ξ, η)ν (t)k = X kπP C ¡ν,j X ¢ FT ρn (t) − j ρW (t)kL1 + |πP L (ν qj )(t) − ν qjW F T (t)|. (28) ν j j First of all the evolution inside suppliers of tangent vectors is the same as for the part densities, i.e. shifts evolves travelling at velocity vj on the j-th supplier, and no increase of norms occur. Such increase may occur in three cases: • Interaction of a wave with an empty queue; • Interaction of a wave with a non empty queue; • Emptying of a queue. Let us start analyzing what happens when a wave interacts with an empty queue, say the j-th queue. Let us call ξ− , ∆ρ− , η − , respectively ξ+ , ∆ρ+ , η + , the values of jump shift, jump and queue shift before, respectively after, the interaction. The shifted wave − interacts with the queue with a time delay given by δt = vξj−1 and the new shift is given vj by ξ+ = vj δt = vj−1 ξ− . Since the queue is empty before the interaction, we have η− = 0, while η + = (∆ρ− − ∆ρ+ )ξ− . Therefore µ ¶ ½ ¾ vj vj |ξ+ ∆ρ+ | + |η + | = |ξ− | |∆ρ+ | + |∆ρ− − ∆ρ+ | ≤ max , 1 |ξ− ∆ρ− |, (29) vj−1 vj−1 where we used the fact that ∆ρ− and ∆ρ+ have the same sign. Let us now consider the interaction with a nonempty queue. Calling ξ− , ∆ρ− the values of shifts and jumps before the interaction and by η + the produced shift in the queue, we simply get: η + = ξ− ∆ρ− , (30) thus the norm of the tangent vector is constant. We now focus on the case of the j-th queue emptying at some time t̄ ∈ [n∆tj,ν , (n+1)∆tj,ν [, i.e the queue buffer has positive values in a left neighborhood of t̄ and vanishes in a right neighborhood of t̄. Then the WFT solution will generate a wave on the j-th supplier starting at time t̄, while the projected numerical solution will generate two waves at times 11 n∆tj,ν and (n + 1)∆tj,ν . More precisely at time (n + 1)∆tj,ν , in a right neighborhood of the left endpoint aj , we have: ½ n fj−1,out x ∈ [aj , aj + vj (t̄ − n∆tj,ν )[ j WFT ρν ((n + 1)∆tj,ν , x) = , µj x ∈ [aj + vj (t̄ − n∆tj,ν ), aj + 2 vj ∆tj,ν ] πP C ¡ν,j n+1 ρ ¢ ½ (x) = ρ∗ µj x ∈ [aj , aj + vj ∆tj,ν [ , x ∈ [aj + ∆tj,ν , aj + 2 vj ∆tj,ν ] n n where ρ∗ is a value in the interval [fj−1,out , µj ] (notice that fj−1,out < µj since we assumed that the queue is emptying). This situation can be detected by tangent vectors if we n n , ρ∗ ) and β = consider the discontinuity (fj−1,out , µj ) as splitted in two parts: α = (fj−1,out (ρ∗ , µj ) considering the shifts ξα = −vj (t̄ − n∆tj,ν ) < 0 and ξβ = vj ((n + 1)∆tj,ν − t̄) > 0. These two shifts are caused because of the approximation and the relative norm of tangent vectors are estimated as follows: n |ξα ∆ρα | = vj (t̄ − n∆tj,ν ) (ρ∗ − fj−1,out ), |ξβ ∆ρβ | = vj ((n + 1)∆tj,ν − t̄) (µj − ρ∗ ), therefore we can write: n |ξα ∆ρα | + |ξβ ∆ρβ | ≤ ∆tj,ν (µj − fj−1,out ). (31) Finally we have: Lemma 5 Assume (23) holds true and consider the tangent vector (ξ, η)ν defined by (26), (27). Then it holds:   X X k(ξ, η)ν (t)k ≤ 2−ν K  T V (ρj,0 ) + max{µj−1 , µj } , j j≥2 where T V (·) indicates the total variation and   µ ¶ Y ½ ¾ N ∆xj,0  vj K = max max ,1 . j vj vj−1 (32) j=2 Proof. The norm of the tangent vector at initial time is zero. Then it may increase due to queues emptying. For the j-th queue, such increase, estimated by (31), is bounded by ∆tj,ν times the strengths of waves interacting with the queue, which in turn is bounded F T . On the other side, from (23), we get (25) thus, from by the total variation of j−1 ρW ν formula (2.10a) and (2.10b) at page 165 of [21], we have: ¯ ¯ ¯ X X X ¯¯ d ν qjW F T ¯ FT T V (j ρW (t)) ≤ T V (ρj,0 ) + (0)¯ . ¯ ν ¯ ¯ dt j j j≥2 12 Moreover, it holds: ¯ ¯ ¯ d ν qW F T ¯ ¯ ¯ j (0)¯ ≤ max{µj−1 , µj }. ¯ ¯ ¯ dt The norm of the tangent vector may further increase due to interactions of waves with queues, but such increase is estimated in (29) and (30). Therefore by (22) we conclude. Since WFT approximations are converging to a solution ¡to the¢ODE-PDE coupled model, as proved in [21], from Lemma 5 we get that also (πP C ν,j ρn , πP L (ν qj )) are converging to a solution. We now want further to estimate the convergence rate. For this purpose, define the convergence error as: XX ¯ ¯ X ν n ν+1 n | qj − qj | Eν (n) = 2−ν ∆xj,0 ¯ν,j ρni − ν+1,j ρni ¯ + j X i ¡ν,j ¢ ¡ν+1,j j X (33) Moreover, we have: X X FT j WFT −(ν+1) Eν (0) = kj ρW (0) − ρ (0)k ≤ 2 ∆xj,0 T V (ρj,0 ). 1 L ν ν+1 (34) kπP C ρ − πP C n ¢ ν+1 n qj |. = n ρ kL1 + |ν qjn − j j j j j ρW F T and j ρW F T in time. We then first estimate the increase of the distance between ν ν+1 ¡ ¢ j ρW F T (0) = π ν+1,j ρ0 can be obtained from j ρW F T (0) = Notice that the initial datum P C ν ν+1 ¡ ¢ πP C ν,j ρ0 by possibly shifting waves with tangent vectors with norms of the order 2−(ν+1) ∆xj,0 , in fact both functions are obtained sampling the same BV function on different subgrids. Then we can control the distance by Lemma 8 of the Appendix. More precisely: X X FT j WFT kj ρW (t) − ρ (t)k + |ν qjW F T (t) − ν+1 qjW F T (t)| ≤ 1 L ν ν+1 j j X FT FT kj ρW (0) − j ρW ν ν+1 (0)kL1 . (35) j Now the convergence error can be written as: FT FT Eν (n) ≤ kj ρW (n∆tν ) − j ρW ν ν+1 (n∆tν )kL1 + + X kπP C j + X kπP C ¡ν,j X |ν qjW F T (n∆tν ) − ν+1 W F T qj (n∆tν )| j X ¢ FT ρn (t) − j ρW (t)kL1 + |πP L (ν qj )(t) − ν qjW F T (t)| ν ¡ν+1,j n ρ ¢ j (t) − j WFT ρν+1 (t)kL1 j + X |πP L (ν+1 qj )(t) − ν+1 W F T qj (t)|. j Last four addenda can be estimated using (28) and Lemma 5. While the first two addenda are estimated using (35) and (34). Finally we get:   X X X Eν (n) ≤ 2−(ν+1) ∆xj,0 T V (ρj,0 )+(2−ν +2−(ν+1) ) K  T V (ρj,0 ) + max{µj−1 , µj } j j 13 j≥2  = 2−ν  1X 2 j   X X 3 ∆xj,0 T V (ρj,0 ) + K  T V (ρj,0 ) + max{µj−1 , µj } , 2 j j≥2 where K is defined in (32). Finally we get the following: Theorem 6 Assume that ρj,0 are BV functions, ρj,0 (x) ≤ µj for every x and (22) holds true. Then the convergence error Eν (n) (defined in (33)) tends to zero uniformly in n with linear convergence rate in ∆xj,ν = 2−ν ∆xj,0 . 5 Numerical tests : Part I We present some simulation results to illustrate the outcome of methods and discussion of previous Sections. To fix notation, let us define: • ordinary method (OM): Upwind scheme for densities, equation (12); Euler scheme for queues, equation (14); • different spatial steps method (DSSM): Upwind scheme for densities, equation (19), with different CFL conditions, equation (20); Euler scheme for queues, equation (14); • different temporal steps method (DTSM): Upwind method for densities, equation (12); modified Euler scheme for queues, equation (21). 5.1 First test - negative queue buffer occupancies In this test we show the importance of flux corrections in case of oscillating queue buffer occupancies (see subsection 3.1). Consider a supply chain with N = 4 suppliers and the following parameters: Lj = Tj = 1, j = 1, .., 4; µ1 = 99, µ2 = 12, µ3 = 10, µ4 = 8. The initial conditions are: ρj,0 ≡ 0, j = 1, ..., 4, qj,0 = 0, j = 2, 3, 4. A total simulationµtime T = 300 ¶ is considered and for the 3πt first arc the inflow profile is given by f1 (t) = 7.5 1 + sin . 20 Figure 3 depicts the evolutions of uncorrected and corrected queue buffer occupancies q2 (t), simulated by OM with ∆x = ∆t = 0.2. It shows that meaningful numerical errors occur if flux corrections are not used. At time t = 270, the relative difference among peaks of corrected and uncorrected queue buffer occupancies is about 29%. In terms of L1 norm the relative error in about 44.4% for t ∈ [270, 300]. An erroneous approximation of q2 (t) has a cascade effect, affecting the evolution of parts density for the supplier 2, but then also q3 (t), q4 (t), ρ3 (t, x) and ρ4 (t, x). In Figure 4 ρ3 (30, x) is depicted with and without correction. Notice that the presence of a negative q2 (t) generated an evolution of ρ3 (30, x) which is completely different from the real dynamic. At x = 0.6, the relative error in ρ3 (30, x) is about 40%. 14 q2 HtL 10 8 6 4 2 0 260 270 280 290 t 300 Figure 3: Queue buffer occupancy q2 : behaviour without flux correction (dashed line) and with flux correction (solid line). ΡH30,xL 10 9 8 7 6 0.52 0.54 0.56 0.58 0.6 x Figure 4: Behaviour of the density on the arc 3 at t = 30 with flux correction (dashed line) and without flux correction (continuous line). 15 5.2 Second test - different space and time meshes We consider a supply chain with N = 4 suppliers. Maximal fluxes, processing times and lengths of each supplier are reported in the following table (see [17]): supplier j 1 2 3 4 µj 25 15 10 15 Tj 1 1 3 1 Lj 1 0.2 0.6 0.2 We assume that all arcs are empty at t = 0, i.e. ρj,0 ≡ 0; also queues at t = 0 are assumed empty: qj,0 = 0, j = 2, 3, 4. The total simulation time is T = 140. The inflow for the first supplier is given by:  18 0 ≤ t ≤ T4 ,  35 t, 18 f1 (t) = (36) 36 − 35 t, T4 < t < T2 ,  T 0, 2 ≤ t ≤ T. In Figure 5 (left), we present the simulation results for queue buffer occupancies, obtained by OM with ∆x = 0.0125, ∆t = 0.5∆x. Different behaviour for queues occur: q4 remains empty as µ4 > µ3 , while q3 is higher than q2 although processing velocities v2 and v3 are the same. We then run the same test for DSSM and DTSM, choosing, respectively, qi HtL 140 120 100 80 60 40 20 0 q2 HtL 8.4 q2 HtL q3 HtL q4 HtL 8.3 q2 HtL, Dt j q2 HtL, Dx j q2 HtL, ordinary 8.2 8.1 20 40 60 80 t 100 120 140 47.8 47.82 47.84 47.86 47.88 47.9 47.92 t Figure 5: Behaviour of queues for OM (left) and comparison with other methods for q2 (right). ∆t = 0.0125 and ∆x = 0.0125. Other parameters for such methods are the same of OM. In Figure 5 (right), the behaviour of the queue buffer occupancy q2 (t) is shown. Different numerical approximations give rise to very similar results. A further analysis on CPU times (measured in seconds and computed by a Pentium 4, CPU 3.20 GHz, RAM 512 Mb) and convergence errors is made to compare methods. The following tables reports the obtained results. ∆x 0.00625 0.0125 0.025 CPU 1.703 0.156 0.062 16 L1 errors 0.001734 0.079465 0.160575 Table 1: CPU times and L1 errors for OM. ∆t 0.00625 0.0125 0.025 ∆x1 0.00625 0.0125 0.025 ∆x2 0.00125 0.0025 0.005 ∆x3 0.00125 0.0025 0.005 ∆x4 0.00125 0.0025 0.005 CPU 1.671 0.140 0.046 L1 errors 0.002911 0.027847 0.132399 Table 2: CPU times and L1 errors using DSSM. ∆x 0.00625 0.0125 0.025 ∆t1 0.0025 0.05 0.1 ∆t2 0.01875 0.00375 0.075 ∆t3 0.0125 0.0025 0.05 ∆t4 0.00625 0.00125 0.025 CPU 0.828 0.140 0.046 L1 errors 0.003435 0.059411 0.127714 Table 3: CPU times and L1 errors using DTSM. We note that all previous tables contain no meaningful differences, for ∆x = 0.025 and ∆x = 0.0125, in terms of CPU times. Differences occur for the value ∆x = 0.00625 for DTSM because of coarser grids used in three suppliers, thus at the price of less precision. As for L1 errors, they are almost the same for OM, DSSM, and DTSM. We conclude that, as expected, the considered methods have almost the same characteristics, as for goodness of approximation versus CPU times and L1 errors. 6 Improvement of CPU times In order to address the simulation of large networks, we aim at improving computational performances. 6.1 Logic pointer approach Recall the Upwind scheme (12). If space and time grid sizes are set according to the j n relation ∆tj = vj ∆xj then j ρn+1 i+1 = ρi . In other words, the density values at time n + 1 are obtained by those at time n by shifting of one position in the space grid. The computational complexity can be highly reduced by using a logic pointer, that allows to avoid the data shift in the vector of densities. For a better comprehension of such strategy, consider a supplier of length L, divided into n segments. At t = 0, we suppose that the supplier is empty. At t = 1, the density value V1 on the first segment is computed and the pointer is set to position 1. At t = 2, densities should be shifted while the new density value for the first segment, V2 , is computed. The physical shift inside the vector is not considered and the pointer is updated to the second position corresponding to V2 (Figure 6), discarding the previous value and considering cyclic order in the vector. This operation is repeated until a time t = n, while at t = n + 1 the new density is allocated in the first segment and the pointer is set back to position 1 and so on. 17 t=2 Density vector V1 V2 0 0 V2 0 V1 Link 0 0 0 Density on the first segment of the link Figure 6: Situation at t = 2. Operations and computational costs for OM and LP method are examined in the following table where # denotes the cardinality. Operations Densities Computational costs ∼ Boundary data P j Operations Densities Logic pointers L # ∆xjj ∼ #J × Boundary data T ∆tj Computational costs ∼ #J ∼ #J ∼ #J × T ∆tj Table 4: Computational costs for OM (left) and LP (right). Notice that the update of densities for each supplier depends on space grid size for OM. Such dependence vanishes using the LP approach, while the cost of logic pointer update is fixed. Therefore, we get a substantial reduction of computational costs as illustrated in Section 7. 6.2 Different indices numeration A concatenated list is used to store all informations about part densities on suppliers, each one identified by a sequential index. Then, informations for each link are obtained by a sequential research requiring a computational time directly proportional to the number N of suppliers. The computational complexity is O(N ) in the worst case, that occurs when the last arc of the list must be found. In order to decrease the computational complexity, a vector, “Heap” whose i−th element identifies the i− th arc by a pointer, is used. In this way, the generic arc i is found by a direct access to the i− th position of the vector, Heap[i]. The numerical complexity of the research operation becomes O (1), independently from the supply chain dimensions. 7 Numerical tests - Part II We make an analysis of CPU times for the approaches described in Section 6. Consider a sequential supply chain with four suppliers and the following parameters: Lj = 18 Tj = 1, j = 1, 2, 3, 4; µ1 = 25, µ2 = 15, µ3 = 10, µ4 = 15. We first use the method OM with different ∆x = ∆t, total simulation time T = 140, empty queues and arcs as initial condition and inflow profile for the first supplier as in (36). CPU times (measured in seconds and evaluated by a Pentium 4, CPU 3.20 GHz, RAM 512 Mb) and convergence orders are then compared with those of LP and DIN approaches. The following table shows the obtained results. Numerical Code OM LP DIN + LP ∆x = 0.00625 CPU = 0.359, γ = 0.982134 CPU = 0.046, γ = 0.982134 CPU = 0.031 γ = 0.982134 ∆x = 0.0125 CPU = 0.109, γ = 0.991122 CPU = 0.031, γ = 0.991122 CPU = 0.015, γ = 0.991122 ∆x = 0.025 CPU = 0.046, γ = 0.946229 CPU = 0.015, γ = 0.946229 CPU = 0.015, γ = 0.946229 Table 5: CPU times and convergence order (γ ) for different numerical approaches. We then considered a supply chain with N = 10 suppliers with the following parameters: Li = Ti = 1, i = 1, ..., 10; µ1 = 100; µ2 = 7; µ3 = 10; µ4 = 8; µ5 = 15; µ6 = 6; µ7 = 18; µ8 = 7; µ9 = 5; µ10 = 9. For OM we used different ∆x = ∆t, total simulation time T = 140, empty queues and arcs and input profile as in (36). CPU times (measured in seconds and evaluated by a Pentium 4, CPU 3.20 GHz, RAM 512 Mb) and convergence orders are again compared with those obtained by LP and DIN approaches. The following table summarizes the results. Numerical Code OM LP DIN + LP ∆x = 0.00625 CPU = 0.953, γ = 0.900810 CPU = 0.140, γ = 0.900810 CPU = 0.093, γ = 0.900810 ∆x = 0.0125 CPU = 0.312, γ = 1.122965 CPU = 0.078, γ = 1.122965 CPU = 0.046, γ = 1.122965 ∆x = 0.025 CPU = 0.094, γ = 1.023330 CPU = 0.046, γ = 1.023330 CPU = 0.031, γ = 1.023330 Table 6: CPU times and convergence order (γ ) for different numerical approaches. Convergence orders and the obtained numerical results for queues and densities are very similar for OM, LP, and DIN + LP. On the contrary, CPU time can decrease of up to 90% for the two last numerical approaches. This is more evident when the number of arcs of the supply chain is very large, as shown from the following simulation case. Consider a supply chain with N = 100 arcs, and the following characteristics: Li = Ti = 1, i = 1, ..., 100; µ1 = 100; µi = 5, i = 2, ..., 100. For OM we chose different ∆x = ∆t, total simulation time T = 140, empty arcs and queues and input profile as in (36). CPU times (measured in seconds and evaluated by a Pentium 4, CPU 2 GHz, RAM 2 GHz) are compared with those obtained by LP and DIN approaches. The following table reports 19 the obtained results. Numerical Code OM LP DIN + LP ∆x = 0.00625 14.406 5.500 0.953 ∆x = 0.0125 5.719 2.906 0.562 ∆x = 0.025 2.515 1.500 0.265 Table 7: CPU times for different numerical approaches. Appendix Let us first focus on the ρj s: a “generalized tangent vector” consists of two components (v, ξ), where v ∈ L1 describes the L1 infinitesimal displacement, while ξ ∈ Rn describes the infinitesimal displacement of discontinuities. A family of piecewise constant functions θ → ρθ , θ ∈ [0, 1], with the same number of jumps say at the points xθ1 < ... < xθM , admits a tangent vector if the following functions are well defined: ρθ+h (x) − ρθ (x) , h→0 h L1 3 v θ (x)= ˙ lim and also the numbers ξβθ = ˙ lim xθ+h − xθβ β , β = 1, ..., M. h Notice that the path θ → ρθ is not differentiable w.r.t. the¤usual differential structure of £ L1 , in fact if ξβθ 6= 0, as h → 0 the ratio ρθ+h (x) − ρθ (x) /h does not converge to any limit in L1 . The L1 -length of the path γ : θ → ρθ is given by: h→0 kγkL1 Z1 ° ° ° ° = °v θ ° M Z1 ¯ ¯¯ ¯ X ¯ θ ¯¯ ¯ θ dθ + ¯ρ (xβ +) − ρ (xβ −)¯ ¯ξβθ ¯ dθ. 1 L 0 (37) β=1 0 According to (37), the L1 -length of a path γ is also equal to the integral of the norm of its tangent vector, defined as follows: k(v, ξ)k = ˙ kvkL1 + M X |∆ρβ | |ξβ | , β=1 where ∆ρβ = ρ(xβ +) − ρ(xβ −) is the jump across the discontinuity xβ . Now, given two piecewise constant functions ρ and ρ0 , call Ω(ρ, ρ0 ) the family of all “differentiable” paths γ : [0, 1] → γ(t) with γ(0) = ρ, γ(1) = ρ0 . To define a Finsler type metric d, we set the distance between ρ and ρ0 to be equal to © ª d(ρ, ρ0 )= ˙ inf kγkL1 , γ ∈ Ω(ρ, ρ0 ) . To define d on all L1 , for given ρ, ρ0 ∈ L1 we set ° ° © d(ρ, ρ0 )= ˙ inf kγkL1 + kρ − ρ̃kL1 + °ρ0 − ρ̃0 °L1 : ª ρ̃, ρ̃0 piecewise constant functions, γ ∈ Ω(ρ, ρ0 ) . It is easy to check that this metric coincides with the usual L1 metric. 20 Remark 7 Since the Finsler type metric d coincides with the usual L1 metric, the reader could think that the whole framework is not so useful. On the contrary, the different differential structure permits to rely on tangent vectors, whose norm can be easily controlled. This would not be possible using the tangent vectors of the usual differential structure of L1 , i.e. having only the v component. Also, while for systems of conservation laws it is possible to find a decreasing functional, this is not the case for networks (see [15]), even for a scalar conservation law. Let us now turn to the supply-chains case. It is easy to see that all paths in L1 connecting piecewise constant functions can be realized using only the ξ component of the tangent vector, see [5]. Therefore, indicating by xβ j the positions of discontinuities, i j = 1, . . . , N , i = 1, . . . , Mj , a tangent vector to a function defined on the network is given by: (ξβ j , ηj ), i where ξβ j is the shift of the discontinuity xβ j , while ηj is the shift of the queue buffer i i occupancy qj . The norm of a tangent vector is given by: X X k(ξβ j , ηj )k = |ξβ j ||∆ρβ j | + |ηj |. i j,i i i j Again, to control the distance among solutions it is enough to control the evolution of norms of tangent vectors. Finally, we have: Lemma 8 The norm of tangent vectors are decreasing along wave front tracking approximations. References [1] D. Armbruster, P. Degond, C. Ringhofer, A model for the dynamics of large queueing networks and supply chains, SIAM Journal on Applied Mathematics, 66 (3), pp. 896 - 920, 2006. [2] D. Armbruster, P. Degond, C. Ringhofer, Kinetic and fluid models for supply chains supporting policy attributes, Transportation Theory Statist. Phys., 2006b. [3] D. Armbruster, D. Marthaler, C. Ringhofer, Kinetic and fluid model hierarchies for supply chains, SIAM J. on Multiscale Modeling, 2 (1), pp. 43 - 61, 2004. [4] A. Bressan, Hyperbolic Systems of Conservation Laws - The One - dimensional Cauchy Problem, Oxford Univ. Press, 2000. [5] A. Bressan, G. Crasta, B. Piccoli, Well-Posedness of the Cauchy Problem for n × n Systems of Conservation Laws, Memoirs of the American Mathematical Society, vol. 146, n. 694 (2000). [6] G.Bretti, R. Natalini, B. Piccoli, Numerical approximations of a traffic flow model on networks, Networks and Heterogeneous Media 1(1), pp. 57-84, 2006. 21 [7] G. Bretti, C. D’Apice, R. Manzo, B. Piccoli, A continuum - discrete model for supply chains dynamics, Networks and Heterogeneous Media (NHM), 2 (4), pp. 661 - 694, 2007. [8] G. Bretti, B. Piccoli, A tracking algorithm for car paths on road networks, SIAM J. on Appl. Dyn. Syst. 7, pp. 510-531, 2008. [9] C. Dafermos, Hyperbolic Conservation Laws in Continuum Physics, Springer - Verlag, 1999. [10] C. Daganzo, A Theory of Supply Chains, Springer Verlag, New York, Berlin, Heidelberg, 2003. [11] C. D’Apice, R. Manzo, A fluid dynamic model for supply chains, Networks and Heterogeneous Media (NHM), 1 (3), pp. 379 - 398, 2006. [12] C. D’Apice, R. Manzo, B. Piccoli, Modelling supply networks with partial differential equations, to appear on Quarterly of Applied Mathematics. [13] D. Helbing, S. Lämmer, T. Seidel, P. Seba, T. Platkowski, Physics, stability and dynamics of supply networks, Physical Review E 70 (2004), 066116. [14] D. Helbing, S. Lämmer, Supply and production networks: from the bullwhip effect to business cycles, in: D. Armbruster, A. S. Mikhailov, and K. Kaneko (eds.) Networks of Interacting Machines: Production Organization in Complex Industrial Systems and Biological Cells, World Scientific, Singapore, 2005, pp. 33 - 66. [15] M. Garavello and B. Piccoli, Traffic flow on networks, Applied Math Series Vol. 1, American Institute of Mathematical Sciences, 2006. [16] S. K. Godunov, A finite difference method for the numerical computation of discontinous solutions of the equations of fluid dynamics, Mat. Sb. 47, 1959, pp. 271 290. [17] S. Göttlich, M. Herty, A. Klar, Network models for supply chains, Communication in Mathematical Sciences, 3(4), pp. 545 - 559, 2005. [18] S. Göttlich, M. Herty, A. Klar, Modelling and optimization of Supply Chains on Complex Networks, Communication in Mathematical Sciences, 4(2), pp. 315 - 330, 2006. [19] S. Göttlich, M. Herty, C. Ringhofer, Optimization of order policies in supply networks, European J. Oper. Res. 202(2), pp. 456-465, 2010. [20] M. Herty, A. Klar, Modeling, Simulation, and Optimization of Traffic Flow Networks, SIAM J. Sci. Comput. 25(3), pp. 1066-1087, 2003 [21] M. Herty, A. Klar, B. Piccoli, Existence of solutions for supply chain models based on partial differential equations, SIAM J. Math. An., 39(1), pp. 160 - 173, 2007. 22 [22] C. Kirchner, M. Herty, S. Göttlich and A. Klar, Optimal control for continuous supply network models, Netw. Heterog. Media 1(4), pp. 675 - 688 ,2006. [23] R. J. Leveque, Finite Volume Methods for Hyperbolic Problems, Cambdrige, 2002. [24] T. Nagatini, D. Helbing, Stability analysis and stabilization strategics for linear supply chains, Physica A: Statistical and Theoretical Physics, 335 (3-4), pp. 644 - 660, 2004. 23