NNew Perspectives on the Erlang-A Queue
Andrew DawSchool of Operations Research and Information EngineeringCornell University257 Rhodes Hall, Ithaca, NY [email protected] PenderSchool of Operations Research and Information EngineeringCornell University228 Rhodes Hall, Ithaca, NY [email protected] 25, 2017
Abstract
The non-stationary Erlang-A queue is a fundamental queueing model that is usedto describe the dynamic behavior of large scale multi-server service systems that mayexperience customer abandonments, such as call centers, hospitals, and urban mobilitysystems. In this paper, we develop novel approximations to all of its transient andsteady state moments, the moment generating function, and the cumulant generatingfunction. We also provide precise bounds for the difference of our approximations andthe true model. More importantly, we show that our approximations have explicitstochastic representations as shifted Poisson random variables . Moreover, we are alsoable to show that our approximations and bounds also hold for non-stationary Erlang-Band Erlang-C queueing models under certain stability conditions. Finally, we performnumerous simulations to support the conclusions of our results.
Keywords:
Multi-Server Queues, Abandonment, Dynamical Systems, Asymptotics, Time-Varying Rates, Moments, Fluid Limits, Erlang-A Queue, Functional Forward Equations,Moment Generating Function, Cumulant Moment Generating Function.
Markov processes are important modeling tools that help researchers describe real-worldphenomena. Thus, it comes as no surprise that the Erlang-A model, which is a Markovianand multi-server queueing model that incorporates customer abandonments, is an important1 a r X i v : . [ m a t h . P R ] D ec odeling tool in a multitude of application settings. Some of the more prominent ap-plications include telecommunications, healthcare, urban mobility and transportation, andmore recently cloud computing. See for example the following work by Mandelbaum et al.[12], Massey [13], Yom-Tov and Mandelbaum [26], Pender and Phung-Duc [24]. Despite itsimportance in many different applications, the Erlang-A queueing model has remained tobe very difficult to analyze and understand. Even the analysis of the moments of Erlang-Aqueue beyond the fourth moment has remained an important topic for additional study.It is well known that the stationary setting of the Erlang-A is much easier to analyze thanits non-stationary counterpart. Some common approaches used to analyze non-stationaryand state dependent queueing models including asymptotic methods such as heavy trafficlimit theory and strong approximations theory, see for example Halfin and Whitt [5], Man-delbaum et al. [11]. Uniform acceleration is extremely useful for approximating the transitionprobabilities and moments such as the mean and variance of Markov processes. Moreover,the strong approximation methods are useful for analyzing the sample path behavior of theMarkov process by showing that the sample paths of properly rescaled queueing processesconverge to deterministic dynamical systems and Gaussian process limits.However, there are two main drawbacks of these asymptotic methods. The first is that themethod is asymptotic as a function of the model parameters and the results really only holdwhen the rates are large and are nearly infinite. Thus, the quality of the approximationsdepends significantly on the size of the model parameters and these asymptotic methodshave been shown to be quite inaccruate for moderate sized model parameter settings, seefor example Massey and Pender [14, 15]. The second main drawback is that the asymptoticmethods do not generate any important insights for the moments or cumulant momentsbeyond order two since the limits are are based on Brownian motion. Since Brownianmotion has symmetry, its cumulants are all zero beyond the second order. Thus, Brownianapproximations are limited in their power to capture asymmetries in higher moments or eventhe dynamics of the moment generating function, cumulant generating function, or Fouriertransform. Moreover, it has been shown recently by Pender [19], Engblom and Pender [2]that the Erlang-A and its variants have non-trivial amounts of skewness and excess kurtosis,which implies that the Erlang-A are not nearly Gaussian for moderate sized queues. Theseresults also demonstrate that it is important to capture the behavior of the Erlang-A modelbeyond its second moment as this information can be used in staffing decisions Massey andPender [16].One common approximation method that is used in the stochastic networks, queueing,and chemical reactions literature is a moment closure approximation . Moment closure ap-proximations are used to approximate the moments of the queueing process with a surrogatedistribution. It is often the case that the set of moment equations for a large number ofqueueing models are not closed, see for example Matis and Feldman [17], Pender [20]. Thus,the closure approximation helps approximate the moments with a closed system using thesurrogate distribution. One such method used by Pender [21], Pender and Ko [22] is to useHermite polynomials for approximating the distribution of the queue length process. In fact,they show that using a quadratic polynomial works quite well. Since the Hermite polyno-mials are orthogonal to the Gaussian distribution, which has support on the entire real line,these Hermite polynomial approximations do not take into account the discreteness of thequeueing process and the fact the queueing process is non-negative. However, they show2hat Hermite polynomials are natural to analyze since they are orthogonal with respect tothe Gaussian distribution and the heavy traffic limits of multi-server queues are Gaussian.In this paper, we perform an in-depth analysis of the moments and the moment generatingfunction of the non-stationary Erlang-A queue. As the Erlang-B and Erlang-C queueingmodels are special cases of the Erlang-A model, we are able to obtain similar results forthose models. Our approach is to use convexity and exploit Jensen’s and the FKG inequalityto obtain bounds on the moments and moment generating function of the Erlang-A queue.What we find even more exciting is that we are able to provide a stochastic representationof our approximations and bounds as a Poisson random variables with a constant shift. Thisshifted Poisson was observed in peer to peer networks by Ferragut and Paganini [3], however,we will show in the sequel, this novel representation will allow us to view our bounds andapproximations in a new way. The main contributions of this work can be summarized as follows: • We provide new approximations for the moments, moment generating function, andcumulant generating function for the nonstationary Erlang-A queue exploting FKGand Jensen’s inequalities. • We derive a novel stochastic intepretation and representation of our approximationsas shifted Poisson random variables or
M/M/ ∞ queues, depending on the context.This sheds new light on the complexity of queues in heavy traffic or critically loadedregimes. • We prove precise error bounds for our approximations and we also prove new upperand lower bounds for the nonstationary Erlang-A queue that become exact in certainparameter settings.
The remainder of this paper is organized as follows. Section 2 introduces the nonstationaryErlang-A queueing model and its importance in stochastic network theory. In Section 3, weprovide approximations for the moments of the Erlang-A system and use these to boundthe true values. In Section 4 we derive approximations for the moment generating functionand cumulant moment generating function of the Erlang-A queue. We again bound the truevalues by these approximations, and we also find a representation for our approximations interms of Poisson random variables or
M/M/ ∞ queues, depending on the context. The Erlang-A queueing model is a fundamental queueing model in the stochastic processesliterature. The work of Mandelbaum et al. [11], shows that the M ( t ) /M/c + M queueing3ystem process Q ≡ { Q ( t ) | t ≥ } is represented by the following stochastic, time changedintegral equation: Q ( t ) = Q (0) + Π (cid:18)(cid:90) t λ ( s ) ds (cid:19) − Π (cid:18)(cid:90) t µ · ( Q ( s ) ∧ c ) ds (cid:19) − Π (cid:18)(cid:90) t θ · ( Q ( s ) − c ) + ds (cid:19) , where Π i ≡ { Π i ( t ) | t ≥ } for i = 1 , , Π transformsit into a non-homogeneous Poisson arrival process with rate λ ( t ) that counts the customerarrivals that occured in the time interval [0,t). A random time change for the Poisson process Π , gives us a departure process that counts the number of serviced customers. We implicitlyassume that the number of servers is c ∈ Z + and that each server works at rate µ . Finally, athe random time change of Π gives us a counting process for the number of customers thatabandon service. We also assume that the abandonment distribution is exponential and therate of abandonments is equal to θ .One of the main reasons that the Erlang-A queueing model has been studied so extensivelyis because several important queueing models are special cases of it. One special case is theinfinite server queue. The infinite server queue can be derived from the Erlang-A queuein two ways. The first way is to set the number of servers to infinity. This precludes anyabandonments since the abandonment rate θ · ( Q ( t ) − c ) + is always equal to zero when thenumber of servers is infinite. The second way to derive the infinite server queue is to setthe service rate µ equal to the abandonment rate θ . When µ = θ , this implies that thesum of the service and abandonment departure processes is equal to a linear function i.e. µ · ( Q ( t ) ∧ c ) + θ · ( Q ( t ) − c ) + = µ · Q ( t ) = θ · Q ( t ). Thus, the Erlang-A queueing modelbecomes an infinite server queue.One of the main and important insights of Halfin and Whitt [5] is that for multi-serverqueueing systems, it is natural to scale up the arrival rate and the number of servers simul-taneously. This scaling known as the Halfin-Whitt scaling and been an important modelingtechnique for modeling call centers in the queueing literature. Since the M ( t ) /M/c + M queueing process is a special case of a single node Markovian service network , we can alsoconstruct an associated, uniformly accelerated queueing process where both the new arrivalrate η · λ ( t ) and the new number of servers η · c are both scaled by the same factor η > Halfin-Whitt scaling for the Erlang-A model, we arrive at the followingsample path representation for the queue length process as Q η ( t ) = Q η (0) + Π (cid:18)(cid:90) t η · λ ( s ) ds (cid:19) − Π (cid:18)(cid:90) t µ · ( Q η ( s ) ∧ η · c ) ds (cid:19) − Π (cid:18)(cid:90) t θ · ( Q η ( s ) − η · c ) + ds (cid:19) = Q η (0) + Π (cid:18)(cid:90) t η · λ ( s ) ds (cid:19) − Π (cid:18)(cid:90) t η · µ · (cid:18) Q η ( s ) η ∧ c (cid:19) ds (cid:19) − Π (cid:32)(cid:90) t η · θ · (cid:18) Q η ( s ) η − c (cid:19) + ds (cid:33) . Halfin-Whitt scaling is defined by simultaneously scaling up the rate of customerdemand (which is the arrival rate) with the number of servers. In the context of call centersthis is scaling up the number of customers and scaling up the number of agents to answerthe phones. In the context of hospitals or healthcare this might be scaling up the numberof patients with the number of beds or nurses. Taking the following limits gives us the fluid models of Mandelbaum et al. [11], i.e.lim η →∞ η Q η ( t ) = q ( t ) a . s . (2.1)where the deterministic process q ( t ), the fluid mean , is governed by the one dimensionalordinary differential equation (ODE) • q ( t ) = λ ( t ) − µ · ( q ( t ) ∧ c ) − θ · ( q ( t ) − c ) + . (2.2)Moreover, if one takes a diffusion limit i.e.lim η →∞ √ η (cid:18) η Q η ( t ) − q ( t ) (cid:19) ⇒ ˜ Q ( t ) (2.3)one gets a diffusion process where the variance of the diffusion is given by the following ODE • Var (cid:104) ˜ Q ( t ) (cid:105) = λ ( t ) + µ · ( q ( t ) ∧ c ) + θ · ( q ( t ) − c ) + − · Var (cid:104) ˜ Q ( t ) (cid:105) · ( µ · { q ( t ) < c } + θ · { q ( t ) ≥ c } ) . (2.4) In addtion to using strong approximations to analyze the queue length process one can alsouse the functional Kolmogorov forward equations as outlined in Massey and Pender [15].The functional forward equations for the Erlang-A model are derived as, • E[ f ( Q ( t ))] ≡ ddt E[ f ( Q ( t )) | Q (0) = q (0)] (2.5)= λ · E [ f ( Q ( t ) + 1) − f ( Q ( t ))] + E [ δ ( Q ( t ) , c ) · ( f ( Q ( t ) − − f ( Q ( t )))] , (2.6)for all appropriate functions f and where δ ( Q ( t ) , c ) = µ · ( Q ( t ) ∧ c ) + θ · ( Q ( t ) − c ) + . Forthe special case where f ( x ) = x , we can derive an ode for the mean queue length process as • E[ Q ( t )] = λ ( t ) − µ · E[( Q ( t ) ∧ c )] − θ · E[( Q ( t ) − c ) + ] . (2.7)The first thing to note is that this equation is not autonomous and one needs to know thedistribution of Q ( t ) a priori in order to compute the expectations on the righthand sideof Equation 2.7. To know the distribution a priori is impossible except in some specialcases like the infinite server setting. However, it is easy to derive simple approximationsfor the mean queue length by making some assumptions on the queue length process. Thisis known as a closure approximation and one common closure approximation method is to5imply take the expectations from outside the function to inside the function. This impliesthat the expectation E [ f ( X )] becomes f ( E [ X ]). This method is known as a mean fieldapproximation in physics and is also known as the deterministic mean approximation ofMassey and Pender [15]. By applying the mean field approximation to Equation 2.7, we canshow that the resulting differential equation is given by the following autonomous ODE • E[ Q f ( t )] = λ ( t ) − µ · (E[ Q f ] ∧ c ) − θ · (E[ Q f ] − c ) + . (2.8)By careful inspection, one can observe that the ode given by the mean field approximationis identical to the fluid limit of 2.2. Moreover, if one simulates the queueing process andcompare it to the mean field limit, one notices an ordering property. For example on the leftof Figure 1, we simulate the Erlang-A queue and compare to the fluid model. We observethat when θ < µ , that the simulated mean is larger than the fluid mean. This is preciselywhat our results predict. Moreover, on the right of Figure 1, we simulate the Erlang-Aqueue and compare to the fluid model when θ > µ and observe that simulated queue lengthis smaller than the fluid limit. Time0 10 20 30 40 Q ueue Leng t h Fluid-MeanSim-Mean
Time0 10 20 30 40 Q ueue Leng t h Fluid-MeanSim-Mean
Figure 1: λ ( t ) = 10 + 2 · sin( t ), µ = 1, Q (0) = 0, c = 10. θ = 0 . θ = 2 (Right).Our goal in this work is to explain the behavior that we observe in Figure 1, which wewill do in the following section. Before concluding our overview of the Erlang-A queueingmodel, we make a brief remark for notational clarity. Remark . Throughout the remainder of this work, we use Q ( t ) to represent the truequeueing process and Q f ( t ) to represent the fluid approximation of it. This fluid approxi-mation is a stochastic process that will be fully described in this work. In fact, in Section 4we use characterize the fluid approximations and use insight from these representations tobound the true queue length from above and below. In this section, we prove when the true moments of the Erlang-A queue are either dominatedor dominates their corresponding fluid limit. We find that the relationship between the ser-6ice rate and the abandonment rate determines whether or not the moment is dominated bythe fluid limit. This section is organized as follows. In Subsection 3.1, we derive inequalitiesfor the true mean of the Erlang-A and its fluid approximation. In Subsection 3.2 we extendthese inequalities to analogous results for the m th moment of the queueing system. Finally,in Subsection 3.3 we provide figures from numerical experiments that demonstrate thesefindings. We begin with analysis of the mean of the Erlang-A queue. Before we proceed, we firstestablish a lemma for comparisons of ordinary differential equations that will be fundamentalto our approach to the results.
Lemma 3.1 (A Comparison Lemma) . Let f : R → R be a continuous function in bothvariables. If we assume that initial value problem • x ( t ) = f ( t, x ( t )) , x (0) = x (3.9) has a unique solution for the time interval [0,T] and • y ( t ) ≤ f ( t, y ( t )) for t ∈ [0 , T ] and y (0) ≤ x (3.10) then x ( t ) ≥ y ( t ) for all t ∈ [0 , T ] .Proof. The the proof of this result is given in Hale and Lunel [4].With this lemma in hand, we can now derive relationships for the fluid limit and thetrue mean. As seen in the proof, these results follow from the application of this differentialequation comparison lemma and the convexity seen in the fluid approximation.
Theorem 3.2.
For the Erlang-A queue, if Q (0) = Q f (0) , then the true mean dominates thefluid limit when θ < µ , the fluid limit dominates the true mean when θ > µ , and the twomeans are equal when θ = µ .Proof. Recall that the true mean satisfies the following differential equation • E[ Q ( t )] = λ ( t ) − µ · E[( Q ∧ c )] − θ · E[( Q − c ) + ]and the fluid limit satisfies the following differential equation • E[ Q f ( t )] = λ ( t ) − µ · (E[ Q f ] ∧ c ) − θ · (E[ Q f ] − c ) + . We can simplify both equations by observing that ( X ∧ c ) + ( X − c ) + = X for any randomvariable X. Thus, we have the following two equations for the true mean and the fluid limit • E[ Q ( t )] = λ ( t ) − θ · E [ Q ] + ( θ − µ ) · E[( Q ∧ c )] • E[ Q f ( t )] = λ ( t ) − θ · E[ Q f ] + ( θ − µ ) · (E[ Q f ] ∧ c ) .
7f we take the difference of the two equations, we obtain the following • E[ Q ( t )] − • E[ Q f ( t )] = λ ( t ) − θ · E [ Q ] + ( θ − µ ) · E[( Q ∧ c )] − λ ( t ) + θ · E[ Q f ] − ( θ − µ ) · (E[ Q f ] ∧ c )= θ · (E[ Q f ] − E [ Q ]) + ( θ − µ ) · (E[( Q ∧ c )] − (E[ Q f ] ∧ c ))Now since the minimum function ( Q ∧ c ) is a concave function, we have that(E[( Q ∧ c )] − (E[ Q ] ∧ c )) ≤ θ < µ • E[ Q ( t )] − • E[ Q f ( t )] ≥ , and for θ > µ • E[ Q ( t )] − • E[ Q f ( t )] ≤ . Finally, for θ = µ , we have that • E[ Q ( t )] − • E[ Q f ( t )] = 0since both differential equations are initialized with the same value and the origin is anequilibrium point for the difference. This completes the proof.As discussed in Section 2, the Erlang-A model is quite versatile in its relation to otherqueueing systems of practical interest. In the two following corollaries, we find that Theo-rem 3.2 can be applied to the Erlang-B and Erlang-C models. Corollary 3.3.
For the Erlang-B queueing model, if Q (0) = Q f (0) , then E[ Q ( t )] ≤ E [ Q f ( t )] for all t ≥ .Proof. This is obvious after noticing that the Erlang-B queue is a limit of the Erlang-Aqueue by letting θ → ∞ . Corollary 3.4.
For the Erlang-C queueing model, if Q (0) = Q f (0) , then E[ Q ( t )] ≥ E [ Q f ( t )] for all t ≥ .Proof. This is obvious after noticing that the Erlang-C queue is an Erlang-A queue with θ = 0. Since µ is assumed to be positive, then we fall into the case where θ < µ and thiscompletes the proof. Remark . Given that we use Jensen’s inequality and the FKG inequality later on inthe paper, we find it important to differentiate them. Here we give an example that setsthe two apart. If we have the following function Q n , then Jensen’s inequality implies that E [ Q n ] ≥ E [ Q ] n . However, FKG implies that E [ Q n ] ≥ E [ Q n − ] · E [ Q ]. We find it interestingthat by iterating the FKG inequality n − .2 Inequalities for the m th Moment
In this subsection we will now extend the previous findings for the mean to higher momentsof the queueing system. Like the result for the mean, this is again built through observationof the convexity in the differential equation of the fluid approximation.
Theorem 3.6.
For the Erlang-A queue and m ∈ Z + , if Q (0) = Q f (0) , then E [ Q m ( t )] ≥ E (cid:2) Q mf ( t ) (cid:3) when θ < µ , E [ Q m ( t )] ≤ E (cid:2) Q mf ( t ) (cid:3) when θ > µ , and E [ Q m ( t )] = E (cid:2) Q mf ( t ) (cid:3) when θ = µ .Proof. We will use proof by induction. For the base case we can apply Theorem 3.2. Now,suppose that the statement holds for j ∈ { , , . . . , m − } . Recall that the m th momentsatisfies • E [ Q m ( t )] = λ ( t )E (cid:34) m (cid:88) j =0 (cid:18) mj (cid:19) Q j ( t ) − Q m ( t ) (cid:35) + E (cid:34)(cid:32) m (cid:88) j =0 (cid:18) mj (cid:19) ( − m − j Q j ( t ) − Q m ( t ) (cid:33) (cid:0) θQ ( t ) − ( θ − µ )( Q ( t ) ∧ c ) (cid:1)(cid:35) = λ ( t ) m − (cid:88) j =0 (cid:18) mj (cid:19) E (cid:2) Q j ( t ) (cid:3) + θ m − (cid:88) j =0 (cid:18) mj (cid:19) ( − m − j E (cid:2) Q j +1 ( t ) (cid:3) + ( θ − µ )E (cid:34)(cid:32) m − (cid:88) j =0 (cid:18) mj (cid:19) ( − m − − j Q j +1 ( t ) (cid:33) ∧ (cid:32) c m − (cid:88) j =0 (cid:18) mj (cid:19) ( − m − − j Q j ( t ) (cid:33)(cid:35) and the approximate autonomous version satisfies • E (cid:2) Q mf ( t ) (cid:3) = λ ( t ) m − (cid:88) j =0 (cid:18) mj (cid:19) E (cid:2) Q jf ( t ) (cid:3) + θ m − (cid:88) j =0 (cid:18) mj (cid:19) ( − m − j E (cid:2) Q j +1 f ( t ) (cid:3) + ( θ − µ ) m − (cid:88) j =0 (cid:18) mj (cid:19) ( − m − − j (cid:0) E (cid:2) Q j +1 ( t ) (cid:3) ∧ E (cid:2) cQ j ( t ) (cid:3)(cid:1) = λ ( t ) m − (cid:88) j =0 (cid:18) mj (cid:19) E (cid:2) Q jf ( t ) (cid:3) + θ m − (cid:88) j =0 (cid:18) mj (cid:19) ( − m − j E (cid:2) Q j +1 f ( t ) (cid:3) + ( θ − µ ) (cid:32) E (cid:34) m − (cid:88) j =0 (cid:18) mj (cid:19) ( − m − − j Q j +1 ( t ) (cid:35) ∧ E (cid:34) c m − (cid:88) j =0 (cid:18) mj (cid:19) ( − m − − j Q j ( t ) (cid:35)(cid:33) • E [ Q m ( t )] − • E (cid:2) Q mf ( t ) (cid:3) = λ ( t ) m − (cid:88) j =0 (cid:18) mj (cid:19) E (cid:2) Q j ( t ) − Q jf ( t ) (cid:3) + θ m − (cid:88) j =0 (cid:18) mj (cid:19) ( − m − j E (cid:2) Q j +1 ( t ) − Q j +1 f ( t ) (cid:3) + ( θ − µ ) (cid:32) E (cid:34)(cid:32) m − (cid:88) j =0 (cid:18) mj (cid:19) ( − m − − j Q j +1 ( t ) (cid:33) ∧ (cid:32) c m − (cid:88) j =0 (cid:18) mj (cid:19) ( − m − − j Q j ( t ) (cid:33)(cid:35) − E (cid:34) m − (cid:88) j =0 (cid:18) mj (cid:19) ( − m − − j Q j +1 ( t ) (cid:35) ∧ E (cid:34) c m − (cid:88) j =0 (cid:18) mj (cid:19) ( − m − − j Q j ( t ) (cid:35)(cid:33) . Because the minimum is a concave function, we have that for any X and Y with real meansE [ X ∧ Y ] ≤ E [ X ] ∧ E [ Y ]. Thus, we have that for θ > µ , • E [ Q m ( t )] − • E (cid:2) Q mf ( t ) (cid:3) ≥ , if θ < µ , • E [ Q m ( t )] − • E (cid:2) Q mf ( t ) (cid:3) ≤ , and if θ = µ , • E [ Q m ( t )] = • E (cid:2) Q mf ( t ) (cid:3) = 0since both differential equations are initialized with the same value, the origin is an equi-librium point for the difference, and all the lower-power terms in the differential equationsfollow this structure, which we know from the inductive hypothesis. Therefore we see thisholds for m , which completes the proof.Again as we have seen for the mean, we can exploit the versatility of the Erlang-A queueto extend these insights to the Erlang-B and Erlang-C models as well. Corollary 3.7.
For the Erlang-B queueing model, if Q (0) = Q f (0) , then E[ Q m ( t )] ≤ E (cid:2) Q mf ( t ) (cid:3) for all t ≥ and m ∈ Z + .Proof. This is obvious after noticing that the Erlang-B queue is a limit of the Erlang-Aqueue by letting θ → ∞ . Corollary 3.8.
For the Erlang-C queueing model, if Q (0) = Q f (0) , then E[ Q m ( t )] ≥ E (cid:2) Q mf ( t ) (cid:3) for all t ≥ and m ∈ Z + .Proof. This is obvious after noticing that the Erlang-C queue is an Erlang-A queue with θ = 0. Since µ is assumed to be positive, then we fall into the case where θ < µ and thiscompletes the proof. 10 .3 Numerical Results In this section we describe numerical results for approximating the moments of the Erlang-Aqueue and examine them relative to our findings. In Figures 2 and 3, we show the firstfour moments of the Erlang-A queue and their respective fluid approximations for cases of θ < µ and θ > µ , respectively. In these plots, we take the arrival rate at time t ≥ λ ( t ) = 10 + 2 sin( t ). We initialize the queue as empty, and we assume that the queueingsystem has c = 10 servers each with exponential service rate µ = 1. We test two differentcases for the abandonment rate: θ = 0 . θ = 2. In these settings, we observe that when θ < µ the fluid approximations are below their corresponding simulated stochastic valuesand that when θ > µ the fluid values are greater than the simulations, and this matches thestatements of Theorems 3 and 3.6.We observe the same relationships in Figures 4 and 5. For these plots we instead set λ ( t ) = 100 + 20 sin( t ) and c = 100 and otherwise use the same values as for Figures 2 and 3.With this increase in the arrival intensity and the number of servers, we see that the gapsbetween the fluid approximations and the simulations are again present, albeit proportionallysmaller. Time0 10 20 30 40 Q ueue Leng t h Fluid-MeanSim-Mean (a) First Moment
Time0 10 20 30 40 Q ueue Leng t h Fluid-2nd-MomentSim-2nd-Moment (b) Second Moment
Time0 10 20 30 40 Q ueue Leng t h Fluid-3rd-MomentSim-3rd-Moment (c) Third Moment
Time0 10 20 30 40 Q ueue Leng t h Fluid-4th-MomentSim-4th-Moment (d) Fourth Moment
Figure 2: λ ( t ) = 10 + 2 · sin( t ), µ = 1, θ = 0 . Q (0) = 0, c = 10.11 ime0 10 20 30 40 Q ueue Leng t h Fluid-MeanSim-Mean (a) First Moment
Time0 10 20 30 40 Q ueue Leng t h Fluid-2nd-MomentSim-2nd-Moment (b) Second Moment
Time0 10 20 30 40 Q ueue Leng t h Fluid-3rd-MomentSim-3rd-Moment (c) Third Moment
Time0 10 20 30 40 Q ueue Leng t h Fluid-4th-MomentSim-4th-Moment (d) Fourth Moment
Figure 3: λ ( t ) = 10 + 2 · sin( t ), µ = 1, θ = 2, Q (0) = 0, c = 10. Time0 10 20 30 40 Q ueue Leng t h Fluid-MeanSim-Mean (a) First Moment
Time0 10 20 30 40 Q ueue Leng t h Fluid-2nd-MomentSim-2nd-Moment (b) Second Moment
Time0 10 20 30 40 Q ueue Leng t h Fluid-3rd-MomentSim-3rd-Moment (c) Third Moment
Time0 10 20 30 40 Q ueue Leng t h Fluid-4th-MomentSim-4th-Moment (d) Fourth Moment
Figure 4: λ ( t ) = 100 + 20 · sin( t ), µ = 1, θ = 0 . Q (0) = 0, c = 100.12 ime0 10 20 30 40 Q ueue Leng t h Fluid-MeanSim-Mean (a) First Moment
Time0 10 20 30 40 Q ueue Leng t h Fluid-2nd-MomentSim-2nd-Moment (b) Second Moment
Time0 10 20 30 40 Q ueue Leng t h Fluid-3rd-MomentSim-3rd-Moment (c) Third Moment
Time0 10 20 30 40 Q ueue Leng t h Fluid-4th-MomentSim-4th-Moment (d) Fourth Moment
Figure 5: λ ( t ) = 100 + 20 · sin( t ), µ = 1, θ = 2, Q (0) = 0, c = 100. Building on what we have found for the moments of the Erlang-A, we can provide similarinequalities for the moment generating function and the cumulant generating function againthrough convexity in the differential equations for the fluid approximations. We provide theseinequalities in Subsections 4.1 and 4.2, respectively. In doing so, we find forms for the fluidapproximations that we can interpret in terms of expectations of other random quantities.Through these recognitions, we characterize the fluid approximations. We describe theserepresentations for systems in steady-state in Subsection 4.3 and for nonstationary systemsin Subsection 4.4. We conclude this section with a variety of demonstrations of these resultsthrough empirical experiments in Subsection 4.5.
Using the functional forward equations Massey and Pender [15], we can show that the mo-ment generating function for the Erlang-A queue satisfies the following partial differential13quation • E (cid:2) e α · Q ( t ) (cid:3) = λ ( t ) · ( e α − · E (cid:2) e α · Q ( t ) (cid:3) + θ · ( e − α − · E (cid:2) Q ( t ) · e α · Q ( t ) (cid:3) (4.11) − ( θ − µ ) · ( e − α − · E (cid:2) ( Q ( t ) ∧ c ) · e α · Q ( t ) (cid:3) (4.12)= λ ( t ) · ( e α − · E (cid:2) e α · Q ( t ) (cid:3) + θ · ( e − α − · ∂M ( t, α ) ∂α (4.13) − ( θ − µ ) · ( e − α − · E (cid:2) ( Q ( t ) ∧ c ) · e α · Q ( t ) (cid:3) . (4.14)Just like the non-autonomous differential equation for the mean in Equation 2.7, wealso cannot directly compute the moment generating function since we do not know thedistribution of the queue length a priori. This is also true for numerical purposes. Unless wecan compute the expectation that includes the minimum function it is impossible to knowthe moment generating function, except in special cases such as the infinite server queueand some cases of the Erlang-B queue. Thus, it is useful to obtain approximations thatare explicit upper or lower bounds for the moment generating function. By using Jensen’sinequality for concave functions, we can approximate the moment generating function withthe following partial differential equation • E (cid:2) e α · Q f ( t ) (cid:3) = λ ( t ) · ( e α − · E (cid:2) e α · Q f ( t ) (cid:3) + θ · ( e − α − · ∂M f ( t, α ) ∂α − ( θ − µ ) · ( e − α − · (cid:0) E (cid:2) Q f ( t ) · e α · Q f ( t ) (cid:3) ∧ E (cid:2) c · e α · Q f ( t ) (cid:3)(cid:1) (4.15) ∂M f ( t, α ) ∂t = λ ( t ) · ( e α − · M f ( t, α ) + θ · ( e − α − · ∂M f ( t, α ) ∂α − ( θ − µ ) · ( e − α − · (cid:18) ∂M f ( t, α ) ∂α ∧ c · M f ( t, α ) (cid:19) . (4.16)The following theorem determines exactly when E (cid:2) e α · Q f ( t ) (cid:3) is a lower or upper bound forthe exact moment generating function of the Erlang-A queue. Theorem 4.1.
For the Erlang-A queue, if Q (0) = Q f (0) , then E (cid:2) e α · Q ( t ) (cid:3) ≥ E (cid:2) e α · Q f ( t ) (cid:3) when θ < µ , E (cid:2) e α · Q ( t ) (cid:3) ≤ E (cid:2) e α · Q f ( t ) (cid:3) when θ > µ , and E (cid:2) e α · Q ( t ) (cid:3) = E (cid:2) e α · Q f ( t ) (cid:3) when θ = µ .Proof. If we take the difference of the two partial differential equations, we obtain the fol-lowing • E (cid:2) e α · Q ( t ) (cid:3) − • E (cid:2) e α · Q f ( t ) (cid:3) = λ ( t ) · ( e α − · E (cid:2) e α · Q ( t ) (cid:3) + θ · ( e − α − · E (cid:2) Q ( t ) · e α · Q ( t ) (cid:3) − ( θ − µ ) · ( e − α − · E (cid:2) ( Q ( t ) ∧ c ) · e α · Q ( t ) (cid:3) − λ ( t ) · ( e α − · E (cid:2) e α · Q f ( t ) (cid:3) − θ · ( e − α − · E (cid:2) Q f ( t ) · e α · Q f ( t ) (cid:3) +( θ − µ ) · ( e − α − · (cid:0) E (cid:2) Q f ( t ) · e α · Q f ( t ) (cid:3) ∧ E (cid:2) c · e α · Q f ( t ) (cid:3)(cid:1) = λ ( t ) · ( e α − · (cid:0) E (cid:2) e α · Q ( t ) (cid:3) − E (cid:2) e α · Q f ( t ) (cid:3)(cid:1) + θ · ( e − α − · (cid:0) E (cid:2) Q ( t ) · e α · Q ( t ) (cid:3) − E (cid:2) Q f ( t ) · e α · Q f ( t ) (cid:3)(cid:1) − ( θ − µ ) · ( e − α − · E (cid:2) ( Q ( t ) ∧ c ) · e α · Q ( t ) (cid:3) +( θ − µ ) · ( e − α − · (cid:0) E (cid:2) Q f ( t ) · e α · Q f ( t ) (cid:3) ∧ E (cid:2) c · e α · Q f ( t ) (cid:3)(cid:1) . (cid:2) ( Q ( t ) ∧ c ) · e α · Q ( t ) (cid:3) = E (cid:2)(cid:0) Q ( t ) · e α · Q ( t ) ∧ c · e α · Q ( t ) (cid:1)(cid:3) ≤ (cid:0) E (cid:2) Q f ( t ) · e α · Q f ( t ) (cid:3) ∧ E (cid:2) c · e α · Q f ( t ) (cid:3)(cid:1) . Thus, we have when θ < µ that • E (cid:2) e α · Q ( t ) (cid:3) − • E (cid:2) e α · Q f ( t ) (cid:3) ≥ , (4.17)when θ > µ • E (cid:2) e α · Q ( t ) (cid:3) − • E (cid:2) e α · Q f ( t ) (cid:3) ≤ , (4.18)and finally when θ = µ , • E (cid:2) e α · Q ( t ) (cid:3) − • E (cid:2) e α · Q f ( t ) (cid:3) = 0 (4.19)since they solve the same partial differential equation. This completes our proof.As with the moments, we can observe these relationships occurring in numerical experi-ments. We provide figures demonstrating this in Subsection 4.5. As a consequence of the findings for the moment generating function, we can also providesimilar inequalities for the cumulant moment generating function. Using Equation 4.11, wehave • log (cid:0) E (cid:2) e α · Q ( t ) (cid:3)(cid:1) ≡ ∂∂t log (cid:0) E (cid:2) e α · Q ( t ) (cid:3)(cid:1) = • E (cid:2) e α · Q ( t ) (cid:3) E [ e α · Q ( t ) ] (4.20)= λ ( t ) · ( e α −
1) + θ · ( e − α − · E (cid:2) Q ( t ) · e α · Q ( t ) (cid:3) E [ e α · Q ( t ) ] − ( θ − µ ) · ( e − α − · E (cid:2) ( Q ( t ) ∧ c ) · e α · Q ( t ) (cid:3) E [ e α · Q ( t ) ] (4.21)= λ ( t ) · ( e α −
1) + θ · ( e − α − · ∂G ( t, α ) ∂α − ( θ − µ ) · ( e − α − · E (cid:2) ( Q ( t ) ∧ c ) · e α · Q ( t ) (cid:3) E [ e α · Q ( t ) ] . (4.22)Like for the MGF, we note that we cannot compute the cumulant moment generating functiondirectly without knowing the distribution of the queue length. So, by again applying Jensen’s15nequality, we can describe the fluid approximation as follows. • log (cid:0) E (cid:2) e α · Q f ( t ) (cid:3)(cid:1) = λ ( t ) · ( e α −
1) + θ · ( e − α − · ∂G f ( t, α ) ∂α − ( θ − µ ) · ( e − α − · (cid:32) E (cid:2) Q f ( t ) · e α · Q f ( t ) (cid:3) ∧ E (cid:2) c · e α · Q f ( t ) (cid:3) E [ e α · Q ( t ) ] (cid:33) (4.23) ∂G f ( t, α ) ∂t = λ ( t ) · ( e α −
1) + θ · ( e − α − · ∂G f ( t, α ) ∂α − ( θ − µ ) · ( e − α − · (cid:18) ∂G ( t, α ) ∂α ∧ c (cid:19) . (4.24)Using this observation and our approach in finding the inequalities for the moment gen-erating function, we find the equivalent inequalities for the cumulant moment generatingfunction in the following corollary. Corollary 4.2.
For the Erlang-A queue, if Q (0) = Q f (0) , then log (cid:0) E (cid:2) e α · Q ( t ) (cid:3)(cid:1) ≥ log (cid:0) E (cid:2) e α · Q f ( t ) (cid:3)(cid:1) when θ < µ , log (cid:0) E (cid:2) e α · Q ( t ) (cid:3)(cid:1) ≤ log (cid:0) E (cid:2) e α · Q f ( t ) (cid:3)(cid:1) when θ > µ , and log (cid:0) E (cid:2) e α · Q ( t ) (cid:3)(cid:1) =log (cid:0) E (cid:2) e α · Q f ( t ) (cid:3)(cid:1) when θ = µ .Proof. The proof follows from the same argument that was given in Theorem 4.1 and thefact that the log function is strictly increasing.
From what we have observed for the moment generating function, we can derive an exactrepresentation for the fluid approximation of the moment generating function in steady-state. We assume a stationary arrival rate λ >
0. We will investigate the stationary fluidapproximation differential equations in a casewise manner based on the relationship of λ and the system’s service parameters. To do so, we begin with a lemma bounding the fluidapproximation of the mean. Lemma 4.3.
Suppose that λ is constant. If λ < cµ , then E [ Q f ( ∞ )] < c . Moreover, if λ ≥ cµ , then E [ Q f ( ∞ )] ≥ c .Proof. We will prove this by contradiction. For the first part, we assume that E [ Q f ( ∞ )] ≥ c .Now by using the differential equation for the mean in steady state, we have that0 = λ − µ · ( E [ Q f ( ∞ )] ∧ c ) − θ · ( E [ Q f ( ∞ )] − c ) + = λ − µ · c − θ ( E [ Q f ( ∞ )] − c ) + . Since we assumed that E [ Q f ( ∞ )] ≥ c , then this yields the following inequality λ ≥ cµ, λ ≥ cµ and E [ Q f ( ∞ )] < c , then by the same differential equation we have that λ = µ · ( E [ Q f ( ∞ )] ∧ c ) + θ · ( E [ Q f ( ∞ )] − c ) + = µ · ( E [ Q f ( ∞ )] ∧ c )= cµ + µ · ( E [ Q f ( ∞ )] − c ) < cµ, which yields another contradiction.We now begin characterizing the fluid approximations with our first case, λ ≥ cµ , in thefollowing proposition. Proposition 4.4. If λ ≥ cµ , then in steady-state we have that ∂M f ( ∞ , α ) ∂α = λ · ( e α −
1) + ( θ − µ ) · (1 − e − α ) · cθ · (1 − e − α ) · M f ( ∞ , α ) (4.25) which yields a solution of M f ( ∞ , α ) = e α · ( θ − µ ) · c + λ · ( eα − θ (4.26) for α ∈ R .Proof. To find the partial differential equation, we use functional cumulant bound for anynon-decreasing function h ( · ) (which can be seen as a form of the FKG inequality),E[ h ( X ) · e α · X ]E[ e α · X ] ≥ E[ h ( X )] . (4.27)In the case that λ ≥ cµ we have that E [ Q f ( t )] ≥ c in steady-state by Lemma 4.3, and so weknow how to evaluate the minimum in the fluid equation. Thus, we have that the derivativeof G f ( α ) = log( M f ( ∞ , α )) with respect to α isd G f ( α )d α = λ ( e α −
1) + c ( θ − µ )(1 − e − α ) θ (1 − e − α ) = λe α θ + c ( θ − µ ) θ (4.28)where here we have used the identity e x = e x − − e − x , which can be observed by multiplying eachside of the equation by 1 − e − x . Because the MGF is equal to 1 when α = 0, we also have that G f (0) = 0. Using this initial condition and integrating left and right sides of Equation 4.28with respect to α , we find that G f ( α ) = λ ( e α −
1) + cα ( θ − µ ) θ and since M f ( ∞ , α ) = e G f ( α ) , we attain the stated result.We can now observe that the fluid approximation is equivalent in distribution to a Poissonrandom variable shifted by γ ≡ c ( θ − µ ) θ , as the moment generation function for the Poissondistribution is e β ( e α − , where β is the rate of arrival and α is the space parameter of theMGF. This gives rise to the following. 17 heorem 4.5. For the Erlang-A queue with λ ≥ cµ and m ∈ Z + , if θ > µ E [( Q f ( ∞ ) − γ ) m ] ≤ E [( Q ( ∞ )) m ] ≤ E [( Q f ( ∞ )) m ] and if θ < µ E [( Q f ( ∞ )) m ] ≤ E [( Q ( ∞ )) m ] ≤ E [( Q f ( ∞ ) − γ ) m ] where γ = c ( θ − µ ) θ .Proof. From Proposition 4.4, we have that the fluid approximation of the MGF in steady-state is M f ( ∞ , α ) = e λ ( eα − cα ( θ − µ ) θ = E (cid:2) e α (Γ+ γ ) (cid:3) where Γ ∼ Pois (cid:0) λθ (cid:1) and γ = c ( θ − µ ) θ . From the uniqueness of MGF’s, we have thatE [( Q f ( ∞ )) m ] = E [(Γ + γ ) m ]for all m ∈ Z + . Now, recall that for an M/M/ ∞ queue with arrival rate λ and service rate θ , the stationary distribution is that of a Poisson random variable with rate parameter λθ .So, we can think of Γ as representing the steady-state distribution of an infinite server queuewith Poisson arrival rate λ and exponential service rate θ .Suppose now that θ > µ . Then, by Theorem 3.6 and our preceding observation, we havethat E [( Q ( ∞ )) m ] ≤ E [(Γ + γ ) m ]. Additionally, by comparing the steady-state infinite serverqueue representation of Γ to Q ( ∞ ), we can further observe that E [( Q ( ∞ )) m ] ≥ E [Γ m ], asfor any state j the service rate in Q ( ∞ ) is no more than the service rate in the same statein the Γ queueing system. Thus we have thatE [( Q f ( ∞ ) − γ ) m ] = E [Γ m ] ≤ E [( Q ( ∞ )) m ] ≤ E [(Γ + γ ) m ] = E [( Q f ( ∞ )) m ]for all m ∈ Z + whenever θ > µ . By symmetric arguments, we can also find that if µ > θ then E [( Q f ( ∞ )) m ] = E [(Γ + γ ) m ] ≤ E [( Q ( ∞ )) m ] ≤ E [Γ m ] = E [( Q f ( ∞ ) − γ ) m ]for all m ∈ Z + , as in this case γ = c ( θ − µ ) θ < Remark . Note that in Theorem 3.6, we require that Q (0) = Q f (0) but in this case wehave not assumed such a condition. This is because the inequalities in Theorem 3.6 hold forall time, and we simply need the relationship to hold in steady-state, which can be seen tooccur regardless of initial conditions.By knowing the fluid form of moment generating function explicitly as a Poisson distri-bution, we can also provide exact expressions for the fluid moments and the fluid cumulantmoments. These are given in the two following corollaries. Corollary 4.7. If λ ≥ cµ , then in steady-state we have that the first n moments have thefollowing steady-state expressions: E[ Q nf ( ∞ )] = n (cid:88) j =0 (cid:18) nj (cid:19) · (cid:18) c ( θ − µ ) θ (cid:19) j · P n − j (cid:18) λθ (cid:19) (4.29) where P m (cid:0) λθ (cid:1) is the m th Touchard polynomial with parameter λθ . roof. This can be seen by direct use of the Poisson form of the fluid MGF. Let Γ ∼ Pois (cid:0) λθ (cid:1) and let γ = c ( θ − µ ) θ . Then,E[ Q nf ( ∞ )] = E[(Γ + γ ) n ]= n (cid:88) j =0 (cid:18) nj (cid:19) · γ j · E (cid:2) Γ n − j (cid:3) = n (cid:88) j =0 (cid:18) nj (cid:19) · γ j · P n − j (cid:18) λθ (cid:19) = n (cid:88) j =0 (cid:18) nj (cid:19) · (cid:18) c ( θ − µ ) θ (cid:19) j · P n − j (cid:18) λθ (cid:19) . Corollary 4.8. If λ ≥ cµ , then in steady-state we have that d G f ( ∞ , α )d α (cid:12)(cid:12)(cid:12) α =0 = λθ + c ( θ − µ ) θ = E[ Q f ( ∞ )] (4.30) and for n ∈ Z + d n G f ( ∞ , α )d n α (cid:12)(cid:12)(cid:12) α =0 = λθ = C ( n ) [ Q f ( ∞ )] (4.31) where C ( n ) [ Q f ( ∞ )] is defined as the n th cumulant moment of Q f ( ∞ ) . We now consider the second case, which is λ < cµe − α . Note that this now also requiresa relationship involving the space parameter of the moment generating function, α . This isless general than the first case, but it allows us to derive Lemma 4.9. Lemma 4.9.
For α ≥ , ∂M f ( ∞ , α ) ∂α < cM f ( ∞ , α ) if and only if λ < cµe − α .Proof. To begin, suppose that ∂M f ( ∞ ,α ) ∂α < cM f ( ∞ , α ) . Using this information in conjunctionwith the steady-state form of the partial differential equation for the fluid MGF given inEquation 4.16, we have that0 = λ ( e α − M f ( ∞ , α ) + θ ( e − α − ∂M f ( ∞ , α ) ∂α − ( θ − µ )( e − α − ∂M f ( ∞ , α ) ∂α which simplifies to ∂M f ( ∞ , α ) ∂α = λµ e α M f ( ∞ , α ) . Using our assumption, we see that λµ e α M f ( ∞ , α ) < cM f ( ∞ , α )19nd this yields that λ < cµe − α , which shows one direction.We now move to showing the opposite direction and instead assume that ∂M f ( ∞ ,α ) ∂α ≥ cM f ( ∞ , α ). In this case, Equation 4.16 is given by0 = λ ( e α − M f ( ∞ , α ) + θ ( e − α − ∂M f ( ∞ , α ) ∂α − c ( θ − µ )( e − α − M f ( ∞ , α )and this simplifies to ∂M f ( ∞ , α ) ∂α = λ ( e α −
1) + c ( θ − µ )(1 − e − α ) θ (1 − e − α ) M f ( ∞ , α ) = λe α + c ( θ − µ ) θ M f ( ∞ , α ) . Again by use of this case’s assumption, we have λe α + c ( θ − µ ) θ M f ( ∞ , α ) ≥ cM f ( ∞ , α )and this now yields λ ≥ e − α ( cθ − c ( θ − µ )) = cµe − α , thus completing the proof.We can now use this lemma to find an explicit form for the fluid approximation of thesteady-state moment generating function when λ < cµe − α . Proposition 4.10.
For α ≥ , if λ < cµe − α , then in steady-state we have that ∂M f ( ∞ , α ) ∂α = λ · e α µ · M f ( ∞ , α ) (4.32) which yields a solution of M f ( ∞ , α ) = e λ · ( eα − µ (4.33) for α ∈ R .Proof. By Lemma 4.9 and our assumption that λ < cµe − α , we know that ∂M f ( ∞ ,α ) ∂α Let λ < cµ and m ∈ Z + . Then, if θ > µ E [Γ mθ ] ≤ E [ Q ( ∞ ) m ] ≤ E (cid:2) Γ mµ (cid:3) , (4.34) and if µ > θ E (cid:2) Γ mµ (cid:3) ≤ E [ Q ( ∞ ) m ] ≤ E [Γ mθ ] (4.35) where Γ x ∼ Pois (cid:0) λx (cid:1) for x > .Proof. In each case, the inequality involving Γ µ ∼ Pois (cid:0) λmu (cid:1) follows directly from Proposi-tion 4.10 and Theorem 3.6 via the observation that the fluid form of the moment generatingfunction is equivalent in distribution to that of Γ µ . Here we are using Proposition 4.10 with α = 0, and by continuity we know this holds for some ball around 0. This validates the useof the derivatives of the steady-state MGF with respect to α evaluated at α = 0 in findingthe moments for the fluid approximation. Thus, we are left to prove the inequalities forΓ θ ∼ Pois (cid:0) λθ (cid:1) .To do so, let’s first note that the stationary distribution of a M/M/ ∞ queue with servicerate θ is equivalent to that of Γ θ . Suppose now that θ > µ . Then, any state of such a M/M/ ∞ queue has a larger rate of departure than the same state in the Erlang-A system.Thus, we have that E [Γ mθ ] ≤ E [ Q ( ∞ ) m ] ≤ E (cid:2) Γ mµ (cid:3) for all m ∈ Z + . By symmetric arguments in the θ < µ case, we complete the proof.As we did for the case when λ ≥ cµ , we can use these findings to give explicit expressionsfor the fluid approximations of the moments and the cumulant moments. Corollary 4.13. If λ < cµ , then in steady-state we have that d G f ( ∞ , α )d α (cid:12)(cid:12)(cid:12) α =0 = λµ = E[ Q f ( ∞ )] (4.36) and for n ∈ Z + , d n G f ( ∞ , α )d n α (cid:12)(cid:12)(cid:12) α =0 = λµ = C ( n ) [ Q f ( ∞ )] (4.37)d n M f ( ∞ , α )d n α (cid:12)(cid:12)(cid:12) α =0 = P n (cid:18) λµ (cid:19) = E [ Q f ( ∞ ) n ] (4.38) where C ( n ) [ Q f ( ∞ )] is defined as the n th cumulant moment of Q f ( ∞ ) and P m (cid:16) λµ (cid:17) is the m th Touchard polynomial with parameter λµ . .4 Characterization of the Nonstationary Moment GeneratingFunction Many scenarios that feature customer abandonments may also feature an arrival process thatis nonstationary. To incorporate this, we now incorporate a point process that can be usedto approximate any periodic mean arrival pattern, as discussed in Eick et al. [1]. Specifically,we define λ ( t ) by a Fourier series: let λ and { ( a k , b k ) , k ∈ Z + } be such that λ ( t ) = λ + ∞ (cid:88) k =1 a k sin( kt ) + b k cos( kt ) . (4.39)We now take λ ( t ) as the rate of arrivals at time t in the Erlang-A model. Under this setting,we derive the following expression for the cumulant moment generating function of the fluidapproximation and its corresponding partial differential equation whenever the arrival rateis sufficiently large. We do so through a series of technical lemmas. First, we bound thefluid mean when the arrival rate and initial value are sufficiently large. Lemma 4.14. Suppose that λ ≡ inf t ≥ λ ( t ) > cµ and that E [ Q f (0)] > c . Then, E [ Q f ( t )] > c for all time t ≥ .Proof. We have seen that E [ Q f ( t )] evolves according to • E [ Q f ( t )] = λ ( t ) − µ (E [ Q f ( t )] ∧ c ) − θ (E [ Q f ( t )] − c ) + at all times t . Now, suppose that ˆ t > (cid:2) Q f (ˆ t ) (cid:3) = c + (cid:15) for some (cid:15) > (cid:15) < λ − cµθ we have that • E (cid:2) Q f (ˆ t ) (cid:3) = λ (ˆ t ) − cµ − θ(cid:15) ≥ λ − cµ − θ(cid:15) > . By the continuity of the fluid mean and the fact that E [ Q f (0)] = q (0) > c , we see thatE [ Q f ( t )] > c for all time t ≥ M/M/ ∞ queue with nonstationary arrival rate λ ( t ), which we will use for comparison later in thissection. Lemma 4.15. Let Q ∞ ( t ) be an infinite server queue with nonstationary Poisson arrival rate λ ( t ) and exponential service rate µ and initial value Q ∞ ( t ) = q . Then, E (cid:2) e αQ ∞ ( t ) (cid:3) = e ( e α − (cid:18) λ µ (1 − e − µt )+ (cid:80) ∞ k =1 ( akµ + bkk ) sin( kt )+( bkµ − akk )(cos( kt ) − e − µt ) µ k (cid:19) (cid:0) e − µt ( e α − 1) + 1 (cid:1) q for all t ≥ and α ∈ R . roof. To start, we have that time derivative of the MGF isdE (cid:2) e αQ ∞ ( t ) (cid:3) d t = λ ( t )( e α − (cid:2) e αQ ∞ ( t ) (cid:3) + µ ( e − α − (cid:2) Q ∞ ( t ) e αQ ∞ ( t ) (cid:3) where λ ( t ) is as defined previously: λ ( t ) = λ + ∞ (cid:88) k =1 a k sin( kt ) + b k cos( kt ) . This differential equation can be view as a partial differential equation when expressed as µ (1 − e − α ) ∂M ( α, t ) ∂α + ∂M ( α, t ) ∂t = λ ( t )( e α − M ( α, t )where M ( α, t ) is the moment generating function at time t and space parameter α . Tosimplify our effort, we instead consider the differential equation for the cumulant MGF,which is G ( α, t ) = log( M ( α, t )). This PDE is µ (1 − e − α ) ∂G ( α, t ) ∂α + ∂G ( α, t ) ∂t = λ ( t )( e α − G ( α, 0) = log (cid:0) E (cid:2) e αQ ∞ (0) (cid:3)(cid:1) = log ( e αq ) = αq . Using the notation that G x = ∂G∂x , we seek to solve the system (cid:40) µ (1 − e − α ) G α + G t = λ ( t )( e α − G ( α, 0) = αq and we do so via the method of characteristics. For this approach we introduce variables thecharacteristic variables r and s and establish the characteristic equations, which are ODE’s,as d α d s ( r, s ) = µ (1 − e − α ) , d t d s ( r, s ) = 1 , d g d s ( r, s ) = λ ( t )( e α − α ( r, 0) = r,t ( r, 0) = 0 ,g ( r, 0) = rq . 23e can first see that the ODE’s for α and t solve to α ( r, s ) = log( e c ( r )+ µs + 1) −→ α ( r, s ) = log (( e r − e µs + 1) t ( r, s ) = s + c ( r ) −→ t ( r, s ) = s and so we can now use these to solve the remaining ODE. After substituting we haved g d s ( r, s ) = λ ( s )( e r − e µs which gives a solution of g ( r, s ) = ( e r − (cid:32) λ µ ( e µs − 1) + ∞ (cid:88) k =1 ( a k µ + b k k ) sin( ks ) e µs + ( b k µ − a k k )(cos( ks ) e µs − µ + k (cid:33) + rq . So, using s = t and r = log ( e − µt ( e α − 1) + 1), we have that G ( α, t ) = g (log (cid:0) e − µt ( e α − 1) + 1 (cid:1) , t )= ( e α − (cid:32) λ µ (1 − e − µt ) + ∞ (cid:88) k =1 ( a k µ + b k k ) sin( kt ) + ( b k µ − a k k )(cos( kt ) − e − µt ) µ + k (cid:33) + log (cid:0) e − µt ( e α − 1) + 1 (cid:1) q and therefore by solving for M ( α, t ) = e G ( α,t ) we attain the stated result.Now that we have established these lemmas we proceed with the analysis of the non-stationary Erlang-A. In the next theorem we give explicit forms for the fluid form of thecumulant MGF and its corresponding partial differential equation. Theorem 4.16. If inf t ≤∞ λ ( t ) ≡ λ > cµ and q (0) > c , then in for all t ≥ we have that ∂G f ( t, α ) ∂t = λ ( t ) · ( e α − 1) + θ · ( e − α − · ∂G f ( t, α ) ∂α − c · ( θ − µ ) · ( e − α − 1) (4.40) which gives a solution of G f ( t, α ) = ( e α − (cid:32) λ θ (1 − e − θt ) + ∞ (cid:88) k =1 ( a k θ + b k k ) sin( kt ) + ( b k θ − a k k )(cos( kt ) − e − θt ) θ + k (cid:33) + c ( θ − µ ) θ α + log(( e α − e − θt + 1) (cid:18) q (0) − c ( θ − µ ) θ (cid:19) (4.41) for all t ≥ and all α ∈ R .Proof. From Equation 4.24, we have that the PDE for the fluid approximation’s cumulantmoment generating function is ∂G f ( t, α ) ∂t = λ ( t )( e α − 1) + θ ( e − α − ∂G f ( t, α ) ∂α − ( θ − µ )( e − α − (cid:18) ∂G ( t, α ) ∂α ∧ c (cid:19) . ∂G f ( t,α ) ∂α = E [ Q f ( t ) e αQf ( t ) ] E [ e αQf ( t ) ] . Using the FKG inequality and our observationfrom Lemma 4.14 that E [ Q f ( t )] > c , we have thatE (cid:2) Q f ( t ) e αQ f ( t ) (cid:3) ≥ E [ Q f ( t )]E (cid:2) e αQ f ( t ) (cid:3) > c E (cid:2) e αQ f ( t ) (cid:3) and so (cid:16) ∂G f ( t,α ) ∂α ∧ c (cid:17) = c . Thus, we have the PDE given in Equation 4.40 and so nowwe seek to find it’s solution. We approach this via the method of characteristics. Because G f (0 , α ) = log(E (cid:2) e αQ f (0) (cid:3) ) = αq (0), we see that we seek to solve the following system (cid:40) θ (1 − e − α ) G ( α ) + G ( t ) = λ ( t )( e α − 1) + c ( θ − µ )(1 − e − α ) G f (0 , α ) = αq (0)where G ( x ) = ∂G f ∂x . Introducing characteristic variables r and s , we have the characteristicODE’s as d α d s ( r, s ) = θ (1 − e − α )d t d s ( r, s ) = 1d g d s ( r, s ) = λ ( t )( e α − 1) + c ( θ − µ )(1 − e − α )with initial conditions α ( r, 0) = r , t ( r, 0) = t , and g ( r, 0) = rq (0). Then, we can solve thefirst two ODE’s to see that α ( r, s ) = log(( e r − e θs + 1) t ( r, s ) = s and so we can use these to solve the remaining equation. Substituting in, we have the ODEas d g d s ( r, s ) = λ ( s ) e θs ( e r − 1) + c ( θ − µ ) e θs ( e r − e θs ( e r − 1) + 1and this now solves to g ( r, s ) = ( e r − (cid:32) λ θ ( e θs − 1) + ∞ (cid:88) k =1 ( a k θ + b k k ) sin( ks ) e θs + ( b k θ − a k k )(cos( ks ) e θs − θ + k (cid:33) + c ( θ − µ ) θ (cid:0) log (cid:0) ( e r − e θs + 1 (cid:1) − r (cid:1) + rq (0) . Now, we can rearrange our solutions to find s = t and r = log(( e α − e − θt + 1). Then, wehave that G f ( t, α ) = g (log(( e α − e − θt + 1) , t )= ( e α − e − θt (cid:32) λ θ ( e θt − 1) + ∞ (cid:88) k =1 ( a k θ + b k k ) sin( kt ) e θt + ( b k θ − a k k )(cos( kt ) e θt − θ + k (cid:33) + c ( θ − µ ) θ (cid:0) α − log(( e α − e − θt + 1) (cid:1) + log(( e α − e − θt + 1) q (0)and this simplifies to the stated result. 25ike the approach in our investigation of the steady-state scenario, we can now observethat the fluid approximation is equivalent in distribution to a infinite server queue shiftedby γ ≡ c ( θ − µ ) θ . This gives rise to the following. Theorem 4.17. For the Erlang-A queue with nonstationary arrival rate λ ( t ) such that λ ≡ inf t ≥ λ ( t ) > cµ and initial value q (0) > c , the fluid approximation of the MGF isequivalent to that of a shifted M/M/ ∞ queue with arrival rate λ ( t ) , service rate θ , initialvalue q (0) − c ( θ − µ ) θ , and linear shift c ( θ − µ ) θ .Proof. Observe from Theorem 4.16 that the fluid MGF for the Erlang-A under these condi-tions is M f ( t, α ) = e G f ( t,α ) = e ( e α − (cid:18) λ θ (1 − e − θt )+ (cid:80) ∞ k =1 ( akθ + bkk ) sin( kt )+( bkθ − akk )(cos( kt ) − e − θt ) θ k (cid:19) + c ( θ − µ ) θ α (cid:0) ( e α − e − θt + 1 (cid:1) q (0) − c ( θ − µ ) θ which is of a form that we can recognize. Comparing it to Lemma 4.15, we can see that Q f is of the form of a shifted M/M/ ∞ queue with arrival rate λ ( t ), service rate θ , initialvalue q (0) − c ( θ − µ ) θ , and linear shift c ( θ − µ ) θ , thus enforcing that the fluid model does start at q (0).This representation of the fluid approximation allows us to now provide upper and lowerbounds for the moments of the Erlang-A system. Corollary 4.18. Let Q ( t ) represent the Erlang-A queue with nonstationary arrival rate λ ( t ) such that λ ≡ inf t ≥ λ ( t ) > cµ and initial value q (0) > c , and let Q f ( t ) represent thecorresponding fluid approximation. Then, if θ > µ E [( Q f ( t ) − γ ) m ] ≤ E [ Q ( t ) m ] ≤ E [ Q f ( t ) m ] and if θ < µ E [ Q f ( t ) m ] ≤ E [ Q ( t ) m ] ≤ E [( Q f ( t ) − γ ) m ] for all time t > and all m ∈ Z + , where γ = c ( θ − µ ) θ .Proof. In each case, the bound involving the fluid approximation of the moment is a directconsequence of Theorem 3.6 and so only the other two bounds remain to be shown. Wenow note that since we have characterized the fluid approximation as a shifted M/M/ ∞ queue, the remaining bounds are from the unshifted version of this system and, by followingthe same arguments as in Theorems 4.5 and 4.12 regarding the rates of departure in thecorresponding states of the Erlang-A queue and the M/M/ ∞ queue, this completes theproof. In this subsection we describe various numerical experiments demonstrating these findings.We first have Figures 6, 7, 8, and 9, which compare simulated value of the moment gener-ating function to their fluid approximations. In the first two figures, the arrival intensity26s λ ( t ) = 5 + sin( t ), the service rate is µ = 1, and the number of servers is c = 5. Theabandonment rates are the differing component of these plots, with θ = 0 . θ = 2 as thetwo respective values. These same comparisons are made in the latter two figures, howeverin this case the arrival rate is instead λ ( t ) = 10 + 2 sin( t ) and the number of servers is c = 10.Through these plots one can observe that the true MGF dominates the fluid approxima-tion when θ < µ and that the fluid dominates the stochastic value when θ > µ . This is ofcourse stated with the understanding that for small values of α or for times near 0 the valuesof the MGF and the approximation are quite close and so with numerical error the surfacesmay overlap. Figure 6: λ ( t ) = 5 + sin( t ), µ = 1, θ = 0 . Q (0) = 0, c = 5.Front view (left) and rear view (right).Figure 7: λ ( t ) = 5 + sin( t ), µ = 1, θ = 2, Q (0) = 0, c = 5.Front view (left) and rear view (right).27igure 8: λ ( t ) = 10 + 2 · sin( t ), µ = 1, θ = 0 . Q (0) = 0, c = 10.Front view (left) and rear view (right).Figure 9: λ ( t ) = 10 + 2 · sin( t ), µ = 1, θ = 2, Q (0) = 0, c = 10.Front view (left) and rear view (right).In Figure 10 we plot the limiting distribution for the steady-state Erlang-A. For theseplots we take λ = 20 and µ = 1, and then vary θ and c . For the three plots on the left wetake the abandonment rate to be θ = 0 . θ = 2. For thetop two plots we set the number of servers as c = 15, in the middle two c = 20, and in thebottom two we make c = 25. We observe that the approximate distribution is quite closewhen λ is not near cµ but the approximation is less accurate when λ = cµ . This findingis consistent with much of the literature that focuses on finding novel approximations forqueueing networks and optimal control of these networks, see for example Hampshire andMassey [6], Hampshire et al. [7, 8], Pender and Ko [22], Niyirora and Pender [18], Qin andPender [25]. We note here that these approximations are not all of the same form: recallthat when λ ≥ cµ the fluid approximation is equivalent in distribution to a shifted Poissonrandom variable with parameter λθ , but when λ < cµ it is equivalent to a Poisson distributionwith parameter λµ . 28 a) θ = 0 . c = 15 (b) θ = 2, c = 15(c) θ = 0 . c = 20 (d) θ = 2, c = 20(e) θ = 0 . c = 25 (f) θ = 2, c = 25 Figure 10: Empirical and Fluid Limiting Distributions for λ = 20 and µ = 1.In Figure 11 we examine the limiting distributions for the single server case. In theseplots we set µ = 1 and then vary the arrival rate and the abandonment rate. On all plotson the left we set θ = 0 . θ = 2. Further, in the top pair we make λ = 0 . λ = 1, and in the bottom pair λ = 1 . 2. As in Figure 10, Figure 11shows that our approximations are quite good. Thus, we are able to capture single serverdynamics as well as large-scale multi-server dynamics even though they are quite different.This is even more useful as our approximations are non-asymptotic and don’t rely on scalingthe number of servers. 29 a) θ = 0 . λ = 0 . θ = 2, λ = 0 . θ = 0 . λ = 1 (d) θ = 2, λ = 1(e) θ = 0 . λ = 1 . θ = 2, λ = 1 . Figure 11: Empirical and Fluid Limiting Distributions for c = 1 and µ = 1.In Figures 12, 13, and 14, we take the arrival rate as λ ( t ) = 6 . t ), the servicerate as µ = 1, and the number of servers as c = 5. Because inf t ≥ λ ( t ) > cµ , we usethe characterization of the fluid approximation as a shifted M/M/ ∞ queue and compare thesimulated system, the fluid approximation, and the unshifted M/M/ ∞ . In the first figure weconsider the mean for θ = 1 . θ = 0 . θ = 1 . θ = 0 . 9, respectively.30igure 12: Queue Mean for λ ( t ) = 6 . t ), µ = 1, Q (0) = 6, c = 5. θ = 1 . θ = 0 . λ ( t ) = 6 . t ), µ = 1, θ = 1 . Q (0) = 6, c = 5.Front view (left) and rear view (right).31igure 14: MGF for λ ( t ) = 6 . t ), µ = 1, θ = 0 . Q (0) = 6, c = 5.Front view (left) and rear view (right). In this paper we have investigated the Erlang-A queueing system through comparison tothe fluid approximations of its moments and moment generating function as well as of itscumulants and cumulant moment generating function. Through recognizing the convexity inthe differential equations describing these approximations, we have found fundamental rela-tionships between the values of these quantities and their fluid counterparts: when the rateof abandonment is less than the rate of service the true value dominates the approximation,when the service rate is larger the approximation dominates the true value, and when therates of abandonment and service are equal, the two are equivalent.In forming these inequalities, we have found explicit representations of the fluid approx-imations through equivalences in distribution with Poisson random variables and infiniteserver queues, in the stationary and non-stationary cases, respectively. These characteri-zations both give insight into the approximations themselves and yield natural inequalitiesthat complement those from the approximations. We have demonstrated the performanceof these bounds through simulations. Through consideration of both these findings and theempirical experiments, we can identify interesting directions of future work.For example, it would be of great interest to gain more explicit insights into the gapbetween the fluid approximations and the true values. This is a non-trivial endeavor, whichstems from the non-differentiablility and non-closure in the differential equations for the trueexpectations. The numerical experiments in this work indicate that the fluid approximationsmay often be quite close but not exact, and additional understanding would be useful inpractice. Moreover, extending our results to more complicated queueing systems where thearrival and service processes follow phase type distributions is of interest given the new workof Pender and Ko [22], Ko and Pender [9, 10].Additionally, it would be even more useful to gain a better understanding of the limitingdistribution of the Erlang-A queue. As we discuss in the paper, the empirical experimentsin Subsection 4.5 indicate that the true limiting distributions closely resemble the shifted32oisson distributions that we have found as characterizations of our fluid approximations. Inparticular, the approximations seem quite close when λ is not near cµ . As a simple extensionof this work, it can be observed that some sort of combination of the approximation when λ < cµ and of the approximation when λ > cµ could make a nice choice for approximation ofthe distribution when λ = cµ . In some sense, it is not surprising that these approximationsare similar to the true limiting distribution, as the Erlang-A appears to be a M/M/ ∞ queuewith service rate µ (the approximation when λ < cµ ), when only considering the states upto c , and it also resembles some sort of shifted M/M/ ∞ queue with service rate θ (whichalso describes the approximation when λ ≥ cµ ) for states c + 1 and beyond. Finally, it wouldbe interesting to extend this to networks of Erlang-A queues like in Pender and Massey [23],however, we would have to keep track of the routing probabilities carefully to keep track ofthe convexity/concavity of the rate functions. References [1] Stephen G Eick, William A Massey, and Ward Whitt. M t /G/ ∞ queues with sinusoidalarrival rates. Management Science , 39(2):241–252, 1993.[2] Stefan Engblom and Jamol Pender. Approximations for the moments of nonstationaryand state dependent birth-death queues. arXiv preprint arXiv:1406.6164 , 2014.[3] Andres Ferragut and Fernando Paganini. Content dynamics in P2P networks fromqueueing and fluid perspectives. In Proceedings of the 24th International TeletrafficCongress , page 11. International Teletraffic Congress, 2012.[4] Jack K Hale and Sjoerd M Verduyn Lunel. Introduction to functional differential equa-tions , volume 99. Springer Science & Business Media, 2013.[5] Shlomo Halfin and Ward Whitt. Heavy-traffic limits for queues with many exponentialservers. Operations Research , 29(3):567–588, 1981.[6] Robert C Hampshire and William A Massey. Dynamic optimization with applications todynamic rate queues. In Risk and Optimization in an Uncertain World , pages 208–247.INFORMS, 2010.[7] Robert C Hampshire, Otis B Jennings, and William A Massey. A time-varying call cen-ter design via Lagrangian mechanics. Probability in the Engineering and InformationalSciences , 23(02):231–259, 2009.[8] Robert C Hampshire, William A Massey, and Qiong Wang. Dynamic pricing to con-trol loss systems with quality of service targets. Probability in the Engineering andInformational Sciences , 23(02):357–383, 2009.[9] Young Myoung Ko and Jamol Pender. Strong approximations for time-varying infinite-server queues with non-renewal arrival and service processes. Stochastic Models , 2016.[10] Young Myoung Ko and Jamol Pender. Diffusion limits for the ( M AP t /P h t / ∞ ) N queueing network. Operations Research Letters , 45(3):248–253, 2017.3311] Avi Mandelbaum, William A Massey, and Martin I Reiman. Strong approximations forMarkovian service networks. Queueing Systems , 30(1):149–201, 1998.[12] Avi Mandelbaum, William A Massey, Martin I Reiman, Alexander Stolyar, and BrianRider. Queue lengths and waiting times for multiserver queues with abandonment andretrials. Telecommunication Systems , 21(2):149–171, 2002.[13] William A Massey. The analysis of queues with time-varying rates for telecommunicationmodels. Telecommunication Systems , 21(2):173–204, 2002.[14] William A Massey and Jamol Pender. Poster: skewness variance approximation fordynamic rate multiserver queues with abandonment. ACM SIGMETRICS PerformanceEvaluation Review , 39(2):74–74, 2011.[15] William A Massey and Jamol Pender. Gaussian skewness approximation for dynamicrate multi-server queues with abandonment. Queueing Systems , 75(2-4):243–277, 2013.[16] William A Massey and Jamol Pender. Performance and provisioning analysis for thedynamic rate Erlang-A queue. 2017.[17] Timothy I Matis and Richard M Feldman. Transient analysis of state-dependent queue-ing networks via cumulant functions. Journal of Applied Probability , 38(4):841–859,2001.[18] Jerome Niyirora and Jamol Pender. Optimal staffing in nonstationary service centerswith constraints. Naval Research Logistics (NRL) , 63(8):615–630, 2016.[19] Jamol Pender. Gram-Charlier expansion for time varying multiserver queues with aban-donment. SIAM Journal on Applied Mathematics , 74(4):1238–1265, 2014.[20] Jamol Pender. Laguerre polynomial expansions for time varying multiserver queueswith abandonment, 2014.[21] Jamol Pender. Sampling the functional Kolmogorov forward equations for nonstationaryqueueing networks. INFORMS Journal on Computing , 29(1):1–17, 2016.[22] Jamol Pender and Young Myoung Ko. Approximations for the queue length distribu-tions of time-varying many-server queues. INFORMS Journal on Computing , 29(4):688–704, 2017.[23] Jamol Pender and William A Massey. Approximating and stabilizing dynamic rateJackson networks with abandonment. Probability in the Engineering and InformationalSciences , 31(1):1–42, 2017.[24] Jamol Pender and Tuan Phung-Duc. A law of large numbers for M/M/c/Delay off-setup queues with nonstationary arrivals. In International Conference on Analyticaland Stochastic Modeling Techniques and Applications , pages 253–268. Springer, 2016.[25] Ziyuan Qin and Jamol Pender. Dynamic control for nonstationary queueing networks.2017. 3426] Galit B Yom-Tov and Avishai Mandelbaum. Erlang-R: A time-varying queue with reen-trant customers, in support of healthcare staffing.