KINETIC MODELS AND INTRINSIC TIMESCALES: SIMULATION COMPARISON FOR A 2ND ORDER QUEUEING MODEL

. Kinetic models of stochastic production ﬂows can be expanded into deterministic moment equations and thus approximated with appropriate closures. A second order model for the product density and the product speed has previously been proposed. A systematic analysis comparing simulations of the partial diﬀerential equations (PDE) with discrete event simulations (DES) is performed. Speciﬁcally, factory production is modeled as an M/M/1 queue where the arrival process is a non-homogeneous Poisson process. Three fundamental scenarios for such a time dependent inﬂux are studied: An instant step up/step down of the arrival rate, an exponential step up/step down and periodic variation of the average arrival rate. It is shown that the second order model in general yields signiﬁcant improvements over the ﬁrst order model. Adding diﬀusion into the PDE further improves the agreement in particular for queues with low utilization. The analysis also points to fundamental open issues regarding kinetic models of time dependent agent based simulations. Memory eﬀects and the possibility of resonance in deterministic models are caused by intrinsic timescales of the PDE that are not present in the original stochastic processes.


Introduction.
1.1. First and second order kinetic models for M/M/1 queueing systems. Aggregating stochastic flows of products through a factory or supply chain involves converting the stochastic individual arrival and exit processes of parts into arrival and departure rates characterizing the long time average stochastic process. As a result, the expected behavior of a queuing network is described by a set of ordinary differential equation (ODE) models known as fluid models. Fluid equation models are the deterministic equations replacing the random variables with their means where Work in Progress (WIP) at a particular production unit becomes the continuous state variable of the associated ODE. The appeal of fluid models is that they are deterministic dynamical systems that are well understood even though some important issues related to the stability of queuing networks and the stability of the associated fluid models remain unresolved for multi class queuing systems ( [5], [11], [12]).
Fluid models do not really behave like a fluid, because they still treat every production process separately and hence model the production flow through discrete steps. For long production lines with many steps, it makes sense to treat the production steps as a continuum variable and in that way obtain a genuine fluid dynamical description that treats the factory as a pipe and the parts flowing through the factory as a fluid. In contrast to a real fluid, the spatial variable does not describe physical space, rather it denotes the degree of completion of the part, that is, how far along the part is in the system. Calling x ∈ [0, 1] the degree of completion (where x = 0 and x = 1 denote recent arrivals and departures, respectively), ρ(x, t) ≥ 0 describes the density of parts at stage x at time t, and W (t) = 1 0 ρ(x, t)dx is the total WIP in the factory. If the fluid moves with a velocity v(x, t), then the flux is described as F (x, t) = ρ(x, t)v(x, t) and the production rate of the factory is given by F (1, t). Assuming that the defective products are sorted out after the factory production process, there are no sources or sinks in a factory, thus WIP satisfies the mass conservation law where µ is the overall mean production rate of the factory. By a standard argument of transport equations [18] this integral conservation law is equivalent to a differential conservation law of the form Because v(x, t) ≥ 0, the fluid moves from left to right. Hence, the boundary condition is set as the influx F (0, t) = λ(t); i.e. the local flux at stage zero is the arrival rate of the parts into the factory. Together with an initial WIP profile ρ(x, 0) = ρ 0 (x) this sets up a well-defined transport equation (hyperbolic) problem. The crucial modeling part for these aggregate transport equations is the associated flux model: • A constant velocity leads to a linear wave equation corresponding to constant time between entering the domain and exiting the domain (called the cycle time). • Local congestion effects lead to a local flux at position x and time t that depends on the density around that position x and are typically used in traffic flows, Equation (3) is known as the Lighthill-Whitham model and reflects the fact that drivers slow down as the density of cars around them increases. v 0 is the free velocity, i.e. the speed of a single car on the road, while R is the maximal density that a steady state may assume. • A global flux is used for a model that treats the whole factory as a single queue. As a result, the velocity at position x at time t depends on the global quantity of total WIP. Assuming an M /M /1 queue with arrival rate λ and service rate µ the steady-state of the queue is characterized by relationships between mean cycle time τ = 1 v , expected queue length W and arrival and production rate using Little's law [22]: Equation (4) gives us a useful model for the flux of the transport equation, A continuum model for factory production was developed heuristically along these lines in [3]. Similar approaches have been presented in [15,6,13]. Existence of weak solutions for such integro-differential models were proved in [9,17,26] and extensions to control them have been presented by [8,25]. In [2] an extension of the Lighthill-Whitham model was developed from first principles developing a kinetic theory of the stochastic transport processes in a factory model. The approach follows turbulence or gas-dynamical modeling of transport processes ( [7]). Fundamentally, such an approach is based on a probability density distribution describes the probability to find a particle in an x-interval with a speed in a particular v-interval in a certain time interval. To analyze the resulting kinetic equation, a moment expansion relative to the velocity is used which, together with closure assumptions, reduces the infinite set of moment equations down to a finite set [2]. The evolution equation for the first moment reflects mass conservation and thus is exactly given by Equation (2). The appropriate closure then is to set the flux F (ρ) to the steady state value of the production rate of the M/M/1 queue (Eq. (5)). In [3] it was subsequently shown that simulations of this PDE with an adiabatic influx λ(t), i.e. changing slowly in time, leads to an outflux F (1, t) that agrees well with the expected outflux generated from Discrete Event Simulations (DES).
Closing the moment expansion at second order [2] leads to ∂v(x, t) ∂t Together with the steady-state initial values for ρ(x, 0) and v(x, 0) and boundary conditions imposed on the left boundary, x = 0, equations (7), and (8) are a set of well-posed hyperbolic partial differential equations commonly referred to as the pressureless gas dynamics equations. Eq. (8) is the well known Burgers' equation. It models the advection of the variable v(x, t) that is transported along the characteristics. As a result, once the initial material has left the domain, the solution is completely determined by the values taken on at the left boundary. This turns out to resolve the issue of nonlocal velocities associated with the one-dimensional model: The time it takes a part ρdx to move through the factory is determined at the time that the part joins the end of the queue. If the queue is long, i.e. there is a lot of WIP in front of the new arrival, it will take longer to clear the factory. This corresponds to the behavior of an M /M /1 queue in steady-state: the PASTA (Poisson Arrivals See Time Averages) property suggests that an arriving part of a Poisson arrival (even if the Poisson process is not homogeneous) will find an average queue length W and average velocity v given in Eq. (4).
Hence the second order model is completed with the boundary conditions These assumptions generate a dramatically improved model for a production process: The model captures the nonlinear dependence of the cycle times on WIP, and it correctly relates the steady-state WIP, cycle time, and outflux.
While the second order model has been published for some time, no systematic analysis comparing DES and the PDE simulations for time dependent influxes has been done. This paper fills this gap by studying how much the second order model improves the outflow predictions for non-steady-state (i.e. transient) situations compared to a first order model. Not surprisingly variations in the outflow will be correctly resolved, if they result from variations in the local product density that have little effect on the total WIP and hence on the transport velocity. However, when we study the impact of a single large jump in the influx and the cumulative impact of periodic influx we find, on the one hand, significant improvements over the first order model and, on the other hand, systematic deficiencies. We discuss that those deficiencies relate to the interaction between the timescales of the stochastic process and the timescales inherent in the deterministic dynamical model and are a fundamental feature of the standard kinetic modeling for agent based simulations.
1.2. Time dependent arrival rates. We study three time dependent production flows of the M /M /1 queue. The flows are characterized by a time dependent expectation value for the arrival rate λ(t) at the queue. We develop Discrete Event Simulations in χ [16] with an inhomogeneous Poisson process generating the inflow into the queue. Each of the simulated M /M /1 queuing systems will be compared to the continuum models.
The specific simulation details are as follows: The creation of one event list is called a run and the total number of runs for a given input pattern is called a simulation. In order to derive useable statistics, a histogram is built with the runs in the simulation. The simulated-time interval [0, T ] for the DES is partitioned into equispaced intervals of width ∆t = 0.1. Over each interval [t i , t i+1 ) the number of machine exit events that occurred in this interval is counted and then divided by the total run count of that particular simulation. This histogram is then used as the representative profile for the flow though the queueing system. In each of the following scenarios, the statistics are gathered over a sample size of 40,000 runs per simulation. While in general we list the raw production rate as µ in every formula, in the simulations we always set µ = 1.
1.2.1. Scenario 1: Exponential relaxation. We consider a production scenario, whereby the expectation value of the start rate changes from one constant, λ 1 , for t < T 0 to another constant, λ 2 , for t > T 0 . One way to generate a sequence of arrival events that reflect this scenario is to sample the next inter-arrival time from the appropriate distribution when a new particle arrives at the end of the queue. It is important to note that T 0 is not an event in the DES and hence the start rate will only change at the first event following T 0 . Specifically, if the last arrival event before time T 0 happens at t i , then the next arrival at t i+1 occurs after T 0 . The inter-arrival time between event i and i + 1 is thus sampled from the distribution corresponding to rate λ 1 . All arrival events after t i+1 have an inter-arrival time sampled from the distribution characterized by λ 2 . This setup gives an expected arrival pattern of the form: i.e. the expected arrival rate decays exponentially to the new arrival rate with lim t→∞ λ(t) = λ 2 . Figure (1) provides an illustration of a DES simulation process of the arrival rate λ(t) and the theoretical curves given in Eq. (11) for a typical increasing and decreasing transition. For a factory the exponential relaxation input pattern represents a decision by management to reduce the influx to the new level λ 2 at time T 0 . However the immediate execution of the decision is restricted in practice. This delay could be due to supply contracts that must be fulfilled at the original level or that arrivals are the output of other production facilities which are, for some reason, unable to alter their production schedule until their current WIP is completed. This manifests itself in a slow decay/growth to the new rate level λ 2 as individual production resource streams switch over to this new level in time.
1.2.2. Scenario 2: The step. The step scenario is conceptually identical in setup to the exponential relaxation case with a significant exception: At T 0 the rate λ(t) changes from λ 1 to λ 2 without the wait for an arrival event, i.e. This yields a true step-like behavior in λ(t) as seen in Figure 2(a) for a representative decreasing transition.
In order to produce a valid DES for this scenario we use simulation approaches that are tailored to nonhomogenous Poisson processes such as thinning [19,24].
In a factory, a step type input pattern is generated e.g. when the input stream is a merger of two stochastic production streams and one of them is cut at T 0 . Hence the arrival rate instantaneous changes to the arrival rate of the continuing arrival stream.
1.2.3. Scenario 3: The cycle. The third scenario we consider is a cyclic input. Starting from an initial constant rate λ 1 , after t > T 0 the input pattern λ(t) steadily oscillates with a constant amplitude and frequency between the initial rate λ 1 and a minimum rate λ 2 . Precisely, the arrival rate is given by The cyclic frequency ω is varied to generate DES inputs for various periods T = 2π ω , yet with the same amplitude.
Cyclic behavior is widespread in industry. In many industries, the supply of production resources varies periodically over time. For simplicity, elementary Fourier modes were chosen for the cyclic behavior assuming a spectral decomposition of the input pattern. Figure 2(b) shows the averaged DES and the corresponding λ(t) according to Eq. (13).   14)) for the exponential relaxation scenario and the stepwise transition scenario for up and down transitions between λ = 0.3 and λ = 0.7.

2.1.
The first order PDE model does not create the correct outflux pattern for stepwise changes of the influx. As a baseline we study the two step scenarios for the first order model based on mass conservation Figures 3 compares the DES simulations and the simulations of the partial differential equation (Eq. (14) ) for the stepwise transition scenario and the exponential relaxation scenario for up and down transitions between λ = 0.3 and λ = 0.7. The most glaring feature of the first order model is that it generates an inverse response: large outflux spikes in the opposite direction of the transition. This is a result of the non-local velocity. As mass enters/leaves the system the velocity changes, but it does so at every point x identically and instantaneously. Therefore, the outflux being the product of the density ρ and the velocity v at x = 1 will see an instantaneous increase/decrease based on the increase/decrease of v since the ρ advects with finite speed.
We also notice that the transition downwards is modeled far worse than the transition upwards.
2.2. Second order models are better but not good. For convenience we summarize the full second order model: Comparing the DES in Figure 4 with PDE simulations in Figure 3 shows that the second order model performs much better than the first order model. In particular, the inverse response is gone since the velocity is now also advected and hence strictly local. However there remain several troubling aspects: • In all cases, the outflux of the PDE solution has a delay where it stays on the old outflux λ 1 much longer than the DES. Hence the outflux for the decreasing production rate is too high initially, while in the increasing case it is too low initially. • In the cases of decreasing production, Fig. 4(a) and 4(b), the outflux of the PDE model reaches the new steady states quicker than the DES. However, the resulting underproduction balances the overproduction in the previous item and thus the total outflux of the two simulations is identical. • The stepwise transition scenarios, Figure 4(b) and Figure 4(d) show an outflux pattern that on average looks like the outflux pattern of the exponential relaxation scenario but instead of a continuous curve, we find a piecewise constant outflux, adjusting to the new steady state outflux in steps.
The delayed change of the output for t > T 0 and the piecewise constant output pattern of the stepwise transition scenarios have the same cause: The velocity equation (Burgers' equation) describes advection of the velocity v with the speed determined at the boundary x = 0. Hence all the wip that is in the factory at t = T 0 will travel with the speed v = 1 µ−λ1 and hence the outflux at x = 1 is constant and equal to λ 1 until the charcteristic, emanating from x = 0, t = T 0 reaches the end of the domain at x = 1. For transitions starting at high influxes the travel time through the domain is high and hence the delay is long, for transitions with low influxes, the travel time is short and hence the delay is short.
In the exponential relaxation scenario the influx changes continuously and thus the outflux shows only the intial delay step. For the stepwise transition scenario the influx is a step function and hence the wip changes only when the step has moved through the factory, leading to a repeat of the same piecewise constant and stepwise change of the outflux until the new steady state has been reached.

2.3.
Adding diffusion to the model. The completely hyperbolic second order model, Eq. (15) with a step function influx advects the flux profile forward ( Figure  4(b) and 4(d)) whereas the DES for the same scenario can be approximated by a smooth function. This suggests, that diffusion should be added to the model to smooth out steps and corners.
Diffusion is added for the density variable in the mass conservation PDE leading to By trial and error we determined an optimal diffusion coefficient of D = 0.1 where it seems the benefits of diffusion are maximal. Figure 5 shows what can be achieved: • The steps in the outflux as an initial delay or due to the discontinuous influx are all gone. • For any ramp up scenario the match between DES and PDE simulations is very good. • The ramp down scenarios develop a shock wave in the outflux.

2.4.
Periodically driven hyperbolic PDEs have resonances, periodically changing Poisson processes do not. The PDE and the DES for the exponential and stepwise input pattern dealt entirely with the transient behavior of the queueing system. Taking a cue from linear systems theory we are also interested in the systems response of the M/M/1 queue. Treating the queue as a linear operator (which it clearly is not) we study the outflux (response of the system), its amplitude and phase shift, as a function of periodic inputs, represented by pure harmonic functions.
The input pattern has been described already in Eq(13) and is repeated here:     agrees in time and in value, both oscillations have a mean that agrees with the mean of the influx (0.5). There is a slight difference in the minimum output where the PDE simulation seems to be systematically lower than the DES by about 10%. It is remarkable that even the transient decay to the steady oscillation is matched well.
Following the concepts of systems response, we note that the amplitude of the outflux oscillations in steady state are significantly smaller than the amplitude of the influx oscillations, approximately 0.1 for the outflux, vs. 0.2 for the influx, i.e. the queue, independently of the model, acts as a damper, reducing the variation of the outflux.
Varying the influx amplitudes as well as the frequencies we find that: i) as long as the queue is never close to being overloaded, i.e. λ(t) < µ, the output signal always has one dominant frequency for both, the PDE simulations and the DES and it looks close to a harmonic signal. ii) There is always a constant phase shift determined by the steady state from which the cyclic experiment starts.
Amplitudes however, vary dramatically and the signal to noise ratio of the DES is also strongly affected by the frequencies. For instance, Figure 7 shows the DES outflux for a low frequency (T = 24) and a high frequency forcing (T = 1.2). Clearly the steady state amplitudes of the outflux are much different from each other. Additionally, the high frequency forcing leads to a much more noisy output. Fourier analysis of the outflux confirms that the dominant frequency in this signal is the frequency associated with the influx. However, the signal to noise ratio for this frequency is much smaller than for the low frequency forcing.
We use the utilization u(t) = λ(t) µ as a measure of the overall loading of the factory. Focussing on the amplitudes of the steady state outflux oscillations, we can calculate gain-function-like curves for different mean utilization, showing the amplitude of the outflux scaled by the amplitude of the influx as a function of the forcing frequencies.
We find that for small forcing frequencies (less than 0.2), DES and PDE simulations have almost identical gain functions as illustrated in Figure 8.
For higher forcing frequencies a fundamental difference between the DES and the PDE simulations emerges: A resonance appears in the PDE simulation (see Fig. 8(b) and 8(c)), when the forcing period and the travel time for a wave though the domain match and generate a resonant amplification of the outflux amplitude. No such resonance is present for the periodically varying inhomogeneous Poisson process that represents the influx into the queue. The three figures (Fig. 8) differ by their mean loading: This has two effects: i) for low utilization, the queue has lots of spare production capacity and hence will dampen out any fluctuations in the input to generate an almost constant output corresponding to the mean influx. Thus, although the relative error between the DES and the PDE is large in Fig. 8(c) the actual error is rather small as can be inferred from the scale in the figure. ii) Since the cycle time in steady state is given as τ = 1 µ−λ , and since in our simulations µ = 1, the mean frequency associated with a queueν = 1 τ is simply given asν = 1 −ū. Thus the resonance frequency shifts and depends on the mean utilization.

Numerical details.
In solving the PDEs numerically we implemented a simple upwind finite differencing scheme. This was carried out on the numerically conservative form of the system ∂v(x, t) ∂t + ∂ ∂x Descretization in the spatial parameter x ∈ [0,1] was over an equispaced grid of 101 points giving an internodal length dx = 0.01. The evolution of the solution ρ, v ∈ R 100×1 is thus calculated at each time step k by the following standard equations for the upwind method: where ρv k and 1 2 ρ 2 k are the 101 × 1 flux vectors (arithmetic carried out componentwise), dt k is the time step, and D (1) is the 101 × 101 matrix of first-order finite differences. The time step dt k is calculated from the CFL condition as where the maximum is taken over all possible grid components i = 0, 1, . . . , 100.
The value of C = 0.90 was chosen experimentally as it provided the best balance between solution accuracy and numerical stability.
For the case of diffusion in (19) the scheme (21) is extended in the natural way as Here, D is the experimentally derived diffusion coefficient and D (2) is the 101 × 101 matrix of second order finite differences. Finally, the initial values (e.g. ρ 0 , v 0 , ρv 0 , etc) and left boundary values (e.g. ρ k where we used a trapezoidal integration scheme to calculate the total WIP argument.

3.
Conclusion and open problems. Citing [2], there has been a steady stream of papers that model supply chain dynamics via a hyperbolic continuum model. In addition, many more continuum model descriptions for production lines are based on heuristics following a Lighthill-Whitham-Richards traffic model approach [21,23] all of which implicitly derive justification from this paper. Our study is similar to the celebrated requiem by Daganzo [10] and subsequence resurrection by Aw and Rascle [4] of second order models in highway traffic models. Daganzo noticed that a standard approach to a second order model would lead to unphysical backwards travel for traffic flow and Aw and Rascle showed that by introducing a convective derivative for the pressure term, the model could be fixed. The fundamental reason why this worked is that information along the traffic was traveling with the speed of the product. This is not the case in a queuing theory model since changes in the influx λ(t) are policy changes and hence, due to the nonzero probability of a cycle time of zero, travel with infinite speed through the queueing system as shown in Figure 4(a). Forestier-Coste et al [14] tried modeling M/M/1 queues by fitting DES data to a second order model based on the Aw-Rascle idea with a low oder polynomial expansion of the clearing function. While they succeed in tracking the long term cumulative outflow at a time scale of hundreds of cycle times, they do not show agreement for the PDE and the mean of the DES for short term transients.
As discussed already in [2] the kinetic derivation of the first order model is based on an adiabatic closure and thus is valid only for slowly varying influxes. It is common wisdom that using moment closures that lead to higher order models reduces the modeling error and allows to use influxes that change faster than adiabatically. Unfortunately nobody so far has done the actual comparisons between the discrete event simulation and higher oder kinetic models. The current paper remedies this failure and leads to the following conclusions.
3.1. The second order model improves outflow prediction considerably. We have shown that the second order supply chain model, derived in [2] considerably improves the match between the outflux of the PDE simulations and the DES. Specifically we showed that for several different non-homogeneous Poisson influx scenarios the mean sample path of the DES is much better represented by the second order model than by the first order model. Especially when we add diffusion then the step-up scenarios show an excellent match between PDE outflux and DES outflux.
For cyclically varying inputs the difference between the outflux generated by the PDE model and one generated by the DES reduces to a constant phase shift for an influx that is not near a resonance of the PDE. From a practical point of view in order to use the PDE model as an aggregate model for factory production, the constant phase shift can easily be determined through short time window observations and then used to adjust the PDE outflux accordingly.

3.2.
Kinetic theory and timescales. Some of the problems shown for the second order model of an M/M/1 queueing system are fundamental to the whole concept of kinetic models. Kinetic models in general assume ergodicity: time averages are equal to ensemble averages. This is certainly true for gases where the relaxation time of a specific initial condition is very short and hence all transport equations for gases, electrons, etc typically are in an adiabatically evolving equilibrium unless they are very rarified. This is typically not the case for supply chains or other agent based simulations. Relaxation times in highly loaded queues are very long and hence the initial conditions matter for a very long time.
Kinetic theory for processes with multiple timescales have been discussed before, e.g. using an additional ODE to describe the time evolution of the environment in [20]. However, the experiments that are studied still are characterized by the steady state description of the state variables in a changing environment. Production system modeling here is interested in the time evolution of the expected path of a transient phenomenon.
Nevertheless, the promise of kinetic theory to have a deterministic equation whose solution describes the average behavior of a stochastic system remains desirable. In this case transient aggregate models require the determination of the expected sample path for a variable like e.g. the ouflux of a queue as a function of time, averaged over an ensemble of many experiments of the same type. This goal cannot be reached via the solution of a PDE on a time scale for a hydrodynamic limit, if there is no clear scale separation between the stochastic relaxation time and the time scale on which the non-homogeneous Poisson process evolves.
While there is currently no complete theory to do deal with ensemble averages, in some cases we have solutions: • For instance, for the step scenarios in our examples, the PDE simulations show a delay whereas the DES does not. The difference between the change in outflux for the DES and for the PDE simulations is due to the fact that the PDE has an initial condition that corresponds to a fixed queue length. Thus the time it takes for any changes in the input to manifest themselves in the output is equal to the mean waiting time in the queue. Figure 9 shows good agreement between the queueing time for an M/M/1 queue λ µ(µ−λ) as Figure 9. Delay between the PDE simulation and the DES for the exponential relaxation scenario and the stepwise transition scenario. The asterisks mark simulation results and the curve shows the mean waiting time for an M/M/1 queue as a function of the influx parameter λ.
a function of the influx λ and the simulated delays, measured in the PDE simulations. We have shown previously [1] that averaging the PDE solutions over the correct probability distribution for the initial values, will eliminate the delay problem. We can also use this observation to determine the phase shift of the cyclic input scenario analytically. Note however, that this is a problem that generalizes to any non-homogeneous Poisson process. As the PDE changes the influx at a time t = T 0 , the PDE has a unique and deterministic state given by its solution at the time t = T 0 whereas the DES at that time has a probability distribution of states from which it samples and that define the sample paths for t > T 0 .
• Similarly, if the timescale of the non-homogeneous Poisson process is long enough to allow transients to settle before the influx changes, the match between DES and PDE will be very good. This is illustrated by the transient decay of the oscillating outflux where PDE simulation and DES match very well, once the phase shift due to the delay has been adjusted (Fig. 6). • We have shown that diffusion may be a good modeling tool to introduce some un-modeled stochastic effects. Results are particularly good if the stochastic system is in a state that will relax quickly from an initial condition. Thus ramping up from low utilization is characterized by a high velocity (small cycle time) allowing diffusion to adjust the velocity and thus leading to great matches between DES and PDE simulations. We have explored many different steps up and down and found that diffusion is highly effective, even for steps down, if the utilization is less then 0.5.
However, for steps down from high utilization diffusion leads to shock waves that are not in the DES. The major constraint here is that for high utilization the velocity is very small (high cycle time) and hence diffusion has to be restricted to low values in order to not generate flows that flow back out of the domain. Future work will determine whether asymmetric diffusion may resolve this problem.
Similar ideas shows up for the cyclic forcing: The second order model characterizes a deterministic wave equation that includes a feedback from the outflux at x = 1 to the influx at x = 0 due to the boundary condition for the velocity v(0) = µ 1+W (t) , Eq. (10). Thus it is subject to resonant forcing with all its typical results. At the same time a stochastic transport process over finite lengths does not show any resonance, not even anything similar to stochastic resonances, since typically parts move out of the domain relatively fast and are not subject to enough repeated collisions that would amplify a response.
It is an open issue how to modify a wave equation such that it has a unique wave speed and at the same time is not subject to resonant forcing. Currently we are focussing on the derivation of and the modeling error of a completely new hyperbolic PDE model based on the forward Kolmogorov equations for the M/M/1 queues that resolves some of the issues presented in this paper.