Continuous approximation of $M_t/M_t/1$ distributions with application to production

A single queueing system with time-dependent exponentially distributed arrival processes and exponential machine processes (Kendall notation $M_t/M_t/1$) is analyzed. Modeling the time evolution for the discrete queue-length distribution by a continuous drift-diffusion process a Smoluchowski equation on the half space is derived approximating the forward Kolmogorov equations. The approximate model is analyzed and validated, showing excellent agreement for the probabilities of all queue lengths and for all queuing utilizations, including ones that are very small and some that are significantly larger than one. Having an excellent approximation for the probability of an empty queue generates an approximation of the expected outflow of the queueing system. Comparisons to several well-established approximation from the literature show significant improvements in several numerical examples.


Introduction
We consider the case of a single queue with one server, a FIFO (first-in first-out) service rule and Markovian arrival and departure processes with a time-dependent distribution. This reads in Kendall's notation (M t , M t , 1). Thus the system is specified by random arrival of goods and random service times, where the distribution is dependent on time but the system as a whole is Markovian. Determining the expected outflow of the production is an important performance measure and is based on the idle probability of the queuing system. In the stationary case, the results are well-known but the transient case is difficult to tackle and a system of infinitely many ODEs has to be solved. Specifically the transient behavior of the probability p k (t) having k customers in the queueing system (queue and service) at time t ≥ 0 is for time-independent rates given by a formula consisting of an infinite sum of modified Bessel functions of the first kind; see [11].
Our interest in this problem comes from the concept of production planning and control in manufacturing industries which has been around since the early 20th century and has a long history of research (for reviews see [2,14]). In essence the problem is to determine the input into a production resource to generate a desired output over time. Operationally the problem splits into a forward problem, which estimates the expected output trajectory given a specified input trajectory, and a backward problem which determines the input pattern required to produce a desired output pattern on average over time.
As production starts, availability of parts, machines and workers are all determined by random processes, there is a long tradition to discuss production planning in the context of queueing theory. In particular performance measures such as the average time in system or queue length under long-run steady-state conditions can be derived for e.g. queuing networks of Jackson type [5]. However for many manufacturing processes, notably in the semi-conductor industry, cycle times are long and planning period short relative to the cycle time violating the assumptions of a steady-state approximation and thus leading to the analysis of time-dependent queuing networks.
At the same time the backward problem is an optimization or an optimal control problem: Given a desired output trajectory, find the optimal input function under the constraint that input and output are related via the solution of the forward problem. If the forward problem can be described by an evolution equation (a set of ordinary or partial differential equations) then such problems can be solved using adjoint calculus [15]. In particular, existence and uniqueness as well as controllability of solutions can be proven [6,7,13] in some cases.
This suggests an attempt to model a time dependent queueing system via a continuous description. There have been two different strands of research in this direction in the last 50 years: Newell in a series of papers in the 1960 suggested a diffusion approximation [19,20,21] and postulated a Fokker Planck equation for the cumulative distribution function F (x, t) for the queue length x at time t. He created models for traffic flows though rush hour. A different approach was introduced by [1,9] based on kinetic theory for the probability density f (x, v, t) of finding a particle at position x in the production process considered as a queue and a machine, moving forward with speed v at time t. Boltzmann equations for f and moment equations with different closures lead to transport equations for the density similar to the Lighthill-Whitham-Roberts model [18,22] for traffic flow and to second order moment equations for the velocity of particles moving through the queues [1].
Recently Armbruster et al. [3] performed a systematic analysis comparing simulations of these moment equations with discrete event simulations (DES) for a factory production modeled as an M/M/1 queue with a non-homogeneous Poisson arrival process. They show that while using higher order moment equations improves the model, these transport equations have intrinsic timescales that are not present in the original stochastic processes and thus may lead to fundamentally bad approximations in some cases.
Our research picks up Newell's idea of a Fokker-Planck equation for the queue length. We will derive an approximation for the probability p k (t) of having k ∈ N 0 customers in the system at time t through a continuous variable ρ(x, t) leading to a drift-diffusion equation known as the Smoluchowski equation. We study the relationship of the queue-length probabilities and the solution the Smoluchowski equation, on the half space [0, ∞) with a linear potential. Additionally, we can find the explicit solution to this equation in the transient case with time-independent rates, based on calculations by Smoluchowksi [24,25].
In the case of time-dependent rates, a numerical scheme is provided and compared to the results obtained by the ODE system. The approximate model is compared to the ODE system solution in various examples and compared in the production context to an approximate formula taken from [23].
Overall, having a PDE model that is a good approximation to a queueing system and that even has explicit solutions for some relevant cases allows us to use a wealth of PDE methods to study this and more complicated queueing systems. This paper is organized as follows. In the second section we state basic definitions and results from the queueing theory, which is needed throughout this paper. The third section addresses the formal derivation of the approximate model for an (M t , M t , 1) queueing distribution, followed by explicit solutions to the associated Smoluchowski equation. This section is followed by the numerical treatment of this model, which is used to study the approximate model and compare it to the exact solution in various examples in section 5. In section 6, the connection to the production context is introduced and the approximate model is compared to a well-established approximation from the literature in numerical examples.
2 Definition and Basic Results on (M t , M t , 1) Queues Let λ : R ≥0 → R ≥0 denote the time-dependent arrival rate, and let µ : R ≥0 → R ≥0 denote the time-dependent processing rate. We denote the number of customers in the system at time t ≥ 0 with the random variable L(t) and define p k (t) = P (L(t) = k) as the probability to have k ∈ N 0 customers in the system at time t. Specifically, we study the behavior of the queueing system given by figure 1.
Service Since the queue length follows a birth-death process, we formally obtain the following relation for k ≥ 1: For k = 0, we obtain The situation of time-independent arrival and processing rates is well established and can be found in the standard literature, e.g., [11]. The steady-state distribution in the case ̺ := λ µ < 1 is given by where ̺ is the so-called traffic intensity or utilization of the queueing system. Since the steadystate distribution is a geometric distribution, the expected queue length is given by 3 Continuous Approximation of (M t , M t , 1) Distributions

Derivation of the Approximating Model
Instead of approximating the exact solution p k (t) to Eq. (1), we derive an approximate model that we solve exactly (or numerically). In [19,20,21], the cumulative distribution function F (t, x) of L(t) is approximated by a second-order PDE. We follow this idea, but we model a probability density function (pdf) instead. Specifically, let us assume that ρ : and we want to derive conditions on the functions a, b and ρ 0 such that The coefficient a(t) will describe the mean behavior of the model, i.e., λ(t) > µ(t) implies an increasing size of customers such that the probability density function is expected to move to the right; in the case λ(t) < µ(t), it is exactly the opposite. Since the system is not deterministic, the coefficient b(t) inherits the variance or fluctuations of the model. Assuming lim x→∞ ρ(x, t) = lim x→∞ ρ x (x, t) = 0 we observe conservation of mass: Since we consider a pdf we set We compute the change in the expected number of customers using (1) as Rearranging the terms yields see, e.g., [23]. This implies that the rate of change of the expected number of customers in the system is given by the arrival minus the service rate plus the service rate multiplied with the idle probability p 0 (t).
On the other hand, integration by parts, assuming lim x→∞ xρ( Comparing (7) and (8), leads to to be satisfied at every time t ≥ 0.
and assumption 3.2 corresponds to the expected increase in the number of customers in the system. Assumption 3.2 guarantees the inclusion of the mean transient behavior. Additionally we want the approximate model to be exact in the steady state (̺ < 1) such that eventually occurring errors become damped. This implies that where ρ is the steady-state solution to (4)-(6) and p k is the steady state distribution given in (3). The steady state solution to (4) is given by ρ(x) =Ce a b x . Since ρ(·, t) and ρ are probability density functions, we have , which equals p k if and only if a b = ln(̺) = ln(λ) − ln(µ), leading to . Altogether, is an approximation of the distribution of the queueing model p k (t) where ρ is the solution to (4)-(6) with the coefficients and initial data ρ 0 satisfying
Since we are interested in the transient behavior of the queueing distribution, we assume that we start with a steady-state distribution, see (9), which is in the form of where ln(µ0)−ln(λ0) , as in (10)- (11), and provided λ0 µ0 < 1. We compare the transition to the new steady state given by determined by λ and µ, again provided λ µ < 1. In this case, we can derive the solution to (4)-(6) with (14) explicitly; it reads To calculate the approximated probabilities p A k (t) of p k (t), we have to integrate (15) We additionally discuss an alternative way to calculate the integrals in the following. Using the mean value theorem we reduce the effort to one evaluation at ξ k (t) of ρ. If we can provide a good approximation of ξ k (t), we reduce the numerical costs in the time-dependent case as well. The continuous-model approximation is exact in steady state, with which is equivalent to This motivates the use of which is continuous in ̺(t) and satisfiesξ k (t) ∈ [k, k + 1]. We define the approximatioñ

Numerical Treatment
To compare and validate the continuous approximation (4)-(6) with the result of the system of ordinary differential equations (1)-(2), we need an approximation of p A k (t) and p k (t). The latter is approximated by reducing the infinite ODE system to N ∈ N equations, where a "boundary" condition is set, i.e., we use the mass conservation to close the equations. We have The resulting system is solved with the Matlab ODE solver ode23 1 . To approximate p A k (t), we distinguish two cases. In the case of constant coefficients a, b and starting with a steady-state distribution given by a 0 , b 0 , we can use the analytic formulas (15) and (16) in combination with (18). The evaluation of the standard normal cumulative distribution function is performed with the Matlab function normcdf 2 .
In the case of time-dependent coefficients a, b, we impose in the following a numerical scheme to approximate the solution of (4)- (6). Let {i∆x : i ∈ N 0 } be a discretization of the half space [0, ∞) with fineness ∆x > 0. Since equations (4)-(6) imply the conservation of mass and since equation (5) is a no-flux boundary condition, it is natural to start with a conservative numerical scheme; see [17]. Let ρ j i be the approximation of ρ(i∆x, t j ) for some time t j ≥ 0, and let be the numerical flux function. We define the iteration by In the case λ(t j ) = µ(t j ), which implies a(t j ) = 0, we only observe diffusion with b(t j ) = λ(t j ), and the solution should decrease in this time step. This is the case if we assume the standard stability condition for diffusion equations, see [12], which reads to be satisfied in every iteration. We use the forward difference in time, which implies the resulting scheme to be first-order accurate in time, and from the second discrete and central derivative, we have second-order accuracy in space.

Computational Results
In the following, we numerically examine the continuous approximations p A k (t) andp k (t) of the queue-length distribution p k (t). In the first part, we consider a steady state at time t = 0, which is determined by the rates λ 0 > 0 and µ 0 > 0 satisfying λ0 µ0 < 1, and the system has an abrupt change to the rates λ 1 > 0 and µ 1 > 0. In this case, we derived the analytic expression for p A k (t) andp k (t), see (16) and (15), such that a numerical scheme for the PDE is not necessary, and we avoid errors arising from the scheme.
The second part addresses the use of the numerical scheme (20)- (21), and we analyze the continuous approximation in the transient case, i.e., time-dependent rates.
We introduce the following measures to evaluate the accuracy of the continuous approximations. Let K ∈ N be the number of equations that we want to compare, we define the supremum error We consider a time horizon T = 100 and K = 100 equations in all the cases, and we restrict the ODE system to N = 1000 equations; see (19).

Ramp Up
Generally, we have two types of steps, the ramp up and the ramp down, where we interpret up and down by the value of the traffic intensity ̺. We first consider three types of ramp ups: moderate, strong and very strong ramp ups. A moderate ramp up is determined by λ 0 = 0.5, λ 1 = 0.8 and µ 0 = µ 1 = 1, which corresponds to an increase in the traffic intensity ̺ 0 = 0.5 to ̺ 1 = 0.8. In figure 2 (a) -(b), the values of p k (t), p A k (t) andp k (t) are shown for different time points and k ∈ {0, 3}. Visually, the continuous approximations given by the squares and diamond markers are coincident with the exact model given by the black dots. Some small displacements can be observed at the second and third time points. This is emphasized by the numerical error measures displayed in figure 2(c). As expected from the derivation of the continuous approximation, we have a decay of the error as time evolves since the model is exact in the steady state. The largest errors occur right after the step at t = 0, which is intuitive since it is the time right at the disturbance and also corresponds to the observation that the analytic formulas (15) and (16) are evaluated at singularities for t → 0, which implies errors. Nevertheless, the maximal difference in the exact model is small, of order 10 −3 , as seen in figure 2(c). We relate "small" to the values given in figure 2 (a)-(b).   The last ramp up example that we use is a very strong ramp up from λ 0 = 0.2 to λ 1 = 2 with µ 0 = µ 1 = 1. In this case, the system has no steady state, which is a fundamental assumption in the derivation of the continuous approximation. We again observe a decreasing supremum error in figure 4(d); however, compared to the values given in 4 (a)-(c), we cannot deduce that they are small. Of course, a constant utilization in time greater than one is not a realistic scenario. We analyze the approximation for short times of over utilization (̺(t) ≥ 1) in the time-dependent case later.

Ramp Down
Since we assume a steady state at t = 0, we only consider a moderate and strong ramp down in the following. The moderate ramp down is given by λ 0 = 0.8, λ 1 = 0.5 and µ 0 = µ 1 = 1.
Visually, the approximations are exact as seen in figure 5 (a)-(b), and supported by the error measures in figure 5(c), which are of order 10 −3 . We see a decreasing error in time again; only the approximationp k (t) shows a different behavior initially. One reason for this phenomenon is the evaluations of the standard normal cdf at singular values, which are different in the case of p k (t).

Time-dependent Coefficient: Cyclic
We consider a cyclic time-varying inflow rate, i.e., as has been studied in [28] to approximate the expected outflow of a queueing system with a second-order model of hyperbolic equations. The parameter λ 0 ≥ 0 denotes the lowest and λ 1 ≥ 0 the highest value of the inflow rate, which is periodic with period T Per > 0. We set the production rate as a constant of µ(t) ≡ 1 and, analogously to the step case, study different values for λ 0 and λ 1 . The probabilities p k (t) are again computed with the ODE system, and the approximate values p A k (t) are approximated with the numerical scheme (20)-(21), where we use ∆x = 0.02, x 0 = 0, x 1 = 200 and a time horizon T = 25. The integration of p A k (t) is done with a trapezoidal rule.
We again study a moderate, strong and very strong case and call the case λ 0 = 0.5 and λ 1 = 0.8 the moderate case. In figure 7, we show the results of a simulation for the probability that no and one customer are in the queueing system, respectively. For both periods, T Per = 10 and T Per = 2, visually, the approximations p A 0 (t) and p A 1 (t) are close to the values p 0 (t) and p 1 (t). Table 1 first column shows the maximal difference between the approximation and the ODE result for k = 0, . . . , 100. We observe that the supremum norm of the error increases with a smaller period T Per but remains of the order 10 −3 , which is small compared to the values in figure 7.  In the strong case given by λ 0 = 0.2 and λ 1 = 0.99 corresponding to larger amplitude oscillations, we can find in figure 8 a larger deviation of the approximate model from the ODE system than in the moderate case. Table 1 second column shows that the numerical error measures are of order 10 −2 in this case and are again increasing as T Per decreases. In the very strong cyclic case λ 0 = 0.2 and λ 1 = 2, we have finite time periods in which the utilization is greater or equal to one and we are in the unstable regime. Nevertheless, visually, the approximations p A 0 (t) and p A 1 (t) are close to the ODE system values, which shows the robustness of the approximation with respect to different utilizations; see figure 9. The errors in table 1 third column are again of order 10 −2 .

Application to Production
The previous section addressed the formal derivation and numerical validation of the continuous approximation given by equations (4)- (6). In this section, we consider the (M t , M t , 1) queueing model in a production context and we derive measures to evaluate them. The interpretation of this queueing model in the production is as follows: We assume that parts arrive (from, e.g., orders) randomly with a mean rate λ(t) and are put into a waiting queue, if the production is busy (here, one unit) or into the processor in the case of an idle production unit. The production time is random, with a mean rate µ(t), and the products are fed into the processor from the storage using a FIFO rule. The analysis of the queue length and number of parts in the system has been well-established, and approximations of the expected number of parts are known; see, e.g., [23,27]. In addition to the length of the queue, the outflow of the production is the most important measure in a production context.
The outflow in [t, t + ∆t] is denoted by for some t ≥ 0 and ∆t > 0, which is a P -a.s. finite random variable. We can compute the expected outflow in [t, t + ∆t] as (1) using (1) and (2). In a natural manner, we define the expected outflow at time t as the limit ∆t → 0, i.e., If we compare the latter with (7), we can write the change rate of the expected number of parts in the system as In numerical experiments, it turned out that using (22) in (23) leads to avoidable numerical errors in the continuous approximation case. To calculate the expected number of parts in the system in the case of the continuous approximation, we use the following idea: If we consider the numerical approximation, we use the centered difference to approximate ρ x (k, t) for every k ∈ N.
In [23], a simple approximation of the expected number of parts is derived and is considered in [27] as well. Let us denote by L K (t) the approximation of the expected number of parts at time t ≥ 0; then, in [23], the approximation satisfies the following initial value problem: for some initial value L K 0 ≥ 0, andT is some parameter used to control the transition. As in the examples in [23], we setT = 0 and L K 0 = ̺(0) 1−̺(0) , which is the expected queue length in the steady state determined by ̺(0).
We consider the example in [23], where the inflow rate and the processing rate are given by figure 10. In figure 11 (a), we compare the queue length approximations using the ODE system (1)-(2), the continuous approximation (4)-(6) and the approximation L K (t) (24). Visually, the continuous approximation L A (t) coincides with the ODE system result L(t). Moreover, the maximal absolute difference is 0.0236, whereas the difference between L K (t) and L(t) is 0.5561 in this example. The latter is shown in figure 11 (a). The expected outflow in figure 11 (b) shows similar effects; the error between the expected outflow computed with the ODE system and the continuous approximation is 0.0457 compared to 0.6761 for the error between the expected outflow and the approximation Out A (t).
Finally, we discuss the expected queue length and outflow approximations for a cyclic inflow and processing rate shown in figure 12. There exist over-saturated (̺(t) ≥ 1) and under-  Figure 11: Comparison of the expected queue length and outflow with inflow and processing rate from [23] saturated (̺(t) < 1) time periods, which implies a strong dynamic in time. For the expected queue length and expected outflow, we observe in figure 13 (a)-(b) that the continuous approximation is again powerful and is close to the ODE system result. Specifically, the maximal absolute error for the expected queue length is 0.1293, and in the case of the expected outflow, we have an error of 0.0196. In this example, the simple approximation in [23] fails to capture all the dynamics, resulting in errors of 1.2527 for the expected queue length and 0.8920 in the case of the expected outflow.

Conclusion
We have derived a continuous approximation of the queue length distributions of an (M t , M t , 1) queueing system based on a Fokker-Planck type of partial differential equation under simple assumptions.
It is instructive to compare our model with the heavy traffic model of queuing theory: The probability density given by (4)-(6) and with the choices (10)-(11) corresponds to a reflected Brownian Motion with drift λ − µ and variance 2(µ(t)−λ(t)) ln(µ(t))−ln(λ(t)) . The diffusion limit of singlestation queues from the literature, see, e.g., [5], leads to a queue-length approximation with drift λ − µ and variance λ + min{λ, µ} in the case of exponentially distributed inter-arrival and service times. In the heavy traffic limit ̺ ր 1 both approximations coincide, whereas for ̺ < 1  Figure 13: Comparison of the expected queue length and outflow with cyclic inflow and processing rate our approximation, see (11), leads to a higher variance in the system. We have shown in various numerical examples that our model approximates the original distribution very well and thus provides a new approach to the forward problem of production planning. In addition the Fokker-Planck equation for an initial step at t = 0 can be solved analytically, allowing the study of solutions for the transient behavior of the continuous approximation as well as the original queueing system. We introduced an appropriate numerical discretization scheme for the approximate model to study the fully transient cases and compare the solutions of the PDE to the solution of a truncated ODE-system for the queue-length distributions. The PDE approximation shows excellent agreement with a truncation of the queue-length distributions at a 1000 modes.
Deriving the output from the Fokker-Planck model relates the model production systems. A comparison of the expected outflow of our model, the true expected outflow based on 1000 ODES and and another well established approximation from the literature shows that our model significantly improves on the literature for all cases considered.
There are several open avenues for future work: • Having a PDE model that is a good approximation to a queueing system and that even has explicit solutions for some relevant cases allows us to use a wealth of PDE methods to study this and more complicated queueing systems. Specifically, the production planning problem now becomes an optimal control problem that can be solved via variational methods [15].
• The improvement of the continuous approximation model over heavy traffic models is essentially due to the fit of the approximation to a known stationary distribution of the M/M/1 model. There are more complicated queueing networks that have stationary distributions that would be obvious candidates for the development of such continuum models.
• Other discrete systems, in particular multi-agent systems, often show a mixture of transport and queuing features that are not resolved well in time-dependent and transient cases.
We expect that with this approach we can derive better models for traffic and pedestrian flows [4,8,10]