Transient behavior of fractional queues and related processes

# Transient behavior of fractional queues and related processes

## Abstract

We propose a generalization of the classical M/M/1 queue process. The resulting model is derived by applying fractional derivative operators to a system of difference-differential equations. This generalization includes both non-Markovian and Markovian properties which naturally provide greater flexibility in modeling real queue systems than its classical counterpart. Algorithms to simulate M/M/1 queue process and the related linear birth-death process are provided. Closed-form expressions of the point and interval estimators of the parameters of the proposed fractional stochastic models are also presented. These methods are necessary to make these models usable in practice. The proposed fractional M/M/1 queue model and the statistical methods are illustrated using financial data.

Keywords: Transient analysis, Fractional M/M/1 queue, Mittag–Leffler function, Fractional birth-death process, Parameter estimation, Simulation.

## 1 Introduction

The M/M/1 queue is without a doubt the simplest model for a queue process. It is characterized by arrivals determined by a Poisson process and an independent service time which is negative-exponentially distributed. It is relatively simple and yet the analysis of its transient behavior leads to considerable difficulties. The main source of these difficulties is the presence of a non-absorbing boundary at zero (empty queue). This means that the analysis becomes simpler when we consider models with absorbing boundaries. As a direct result, the state probability of a linear birth-death process, that is the probability that the queue length is at a specific time , has a particularly nice form.

The aim of this paper is to study some related point processes governed by difference-differential equations containing fractional derivative operators. These processes are direct generalizations of the classical M/M/1 queue and the linear birth-death processes. It is well-known that a fractional derivative operator induces a non-Markovian behavior into a system [see 19]. Moreover, parameter estimation and path generation algorithms of these new fractional stochastic models are derived. Note that the proposed fractional point models (with Markovian and non-Markovian properties) are parsimonious which makes them desirable for modeling real-world non-Markovian queueing systems. Observe that fractional point processes driven by fractional difference-differential equations such as the fractional Poisson, the fractional birth, the fractional death, and the fractional birth-death processes have already been gaining attention more recently [see, e.g., 11, 4, 7, 14, 9].

The article is structured as follows. Section 2 presents the explicit construction of the fractional M/M/1 queue starting from the governing equations and a particular subordination relation. The main result derived in this section is the explicit form of the transient state probabilities for each value of the parameter of fractionality. Information regarding the steady-state behavior (stationary behavior) is also highlighted; in particular the fractional process shares the same steady-state behavior as the classical non-fractional case. In Section 3 we develop closed-form estimators (point and interval) for the model parameters in the case of the fractional linear birth-death process. This is preparatory for the similar subsequent analysis applied to the M/M/1 queue (Section 4). The article ends with an application which shows that our constructed estimators perform well in a real-world example.

## 2 Results for a fractional process related to M/M/1 queues

The classical M/M/1 queue process , , that is the queue length in time can be described by the following difference-differential equations governing the state probabilities , :

 ⎧⎪ ⎪⎨⎪ ⎪⎩ddtpk(t)=−(λ+μ)pk(t)+λpk−1(t)+μpk+1(t),k≥1,ddtp0(t)=−λp0(t)+μp1(t),pk(0)=δk,i, (2.1)

where is the initial number of individuals in the queue and is the Kronecker’s delta. In (2.1) and are the entrance and the service rates, respectively.

To arrive at a possible fractional model we consider the Caputo fractional derivative , , with respect to time . If , , where , is the fractional queue with parameter , the generalized difference-differential equations for the state probabilities with arrival rate , service rate and initial customers, read

 ⎧⎪⎨⎪⎩Dαtpαk(t)=−(λ+μ)pαk(t)+λpαk−1(t)+μpαk+1(t),k≥1,Dαtpα0(t)=−λpα0(t)+μpα1(t),pαk(0)=δk,i. (2.2)

First, we will follow Bailey [2, 3] for the derivation of the probabilities , , but adapting the method to take into considerations the presence of the Caputo derivative. The result obtained by Bailey is the so-called classical solution in terms of modified Bessel functions of the first kind. Note however that the derivation of the state probabilities in the classical case can be carried out in several equivalent ways (see for example Champernowne [8], Parthasarathy [15], Abate and Whitt [1]). In the following we will first treat the solution derived by Bailey and then we will use a simpler but lesser known form due to Sharma [17].

We indicate as the probability generating function.

###### Theorem 2.1.

The Laplace transform , , can be written as

 ~Gα(z,s)=sα−1zi+1−(1−z)[a2(s)]i+1[1−a2(s)]−1−λ[z−a1(s)][z−a2(s)],|z|≤1,R(s)>0. (2.3)

where and are the zeros of .

###### Proof.

From (2.2), we can write

 Dαt[Gα(z,t)−pα0(t)]=−(λ+μ)[Gα(z,t)−pα0(t)]+λzGα(z,t). (2.4)

Using the equation on we have

 DαtGα(z,t) =−λGα(z,t)−μGα(z,t)+μpα0(t)+λzGα(z,t)+μz[Gα(z,t)−pα0(t)], (2.5)

and after simplifying, we obtain, for , the Cauchy problem

 {zDαtGα(z,t)=(1−z)[Gα(z,t)(μ−λz)−μpα0(t)],Gα(z,0)=zi. (2.6)

Applying the Laplace transform to (2.6) leads to

 z[sα~Gα(z,s)−sα−1Gα(z,0)]=(1−z)[~Gα(z,s)(μ−λz)−μ~pα0(s)], (2.7)

where . After some simple algebraic calculations we then have

 ~Gα(z,s)=sα−1zi+1−μ(1−z)~pα0(s)zsα−(1−z)(μ−λz),|z|≤1,R(s)>0. (2.8)

As the above function converges in , the zeros of the numerator and the denominator should coincide. Let us indicate the zeros of the numerator as

 a12(s)=sα+λ+μ±[(sα+λ+μ)2−4λμ]1/22λ, (2.9)

with , . Note that

 ⎧⎨⎩a1(s)+a2(s)=(sα+λ+μ)/λ,a1(s)a2(s)=μ/λ,−λ[1−a2(s)][1−a1(s)]=sα. (2.10)

By Rouché theorem [10, Page 168] we have that the only zero in the unit circle is . Therefore it follows that

 sα−1[a2(s)]i+1−μ[1−a2(s)]~pα0(s)=0, (2.11)

which gives

 ~pα0(s)=sα−1[a2(s)]i+1μ[1−a2(s)]. (2.12)

Now, by considering that

 zsα−(1−z)(μ−λz)=−λ[z−a1(s)][z−a2(s)], (2.13)

equation (2.8) can be rewritten as

 ~Gα(z,s)=sα−1zi+1−(1−z)[a2(s)]i+1[1−a2(s)]−1−λ[z−a1(s)][z−a2(s)]. (2.14)

In the following Theorem 2.2 we prove a subordination relation for the fractional queue , , . This is essential for our next results. Before that, let us introduce some facts on the -stable subordinator and its inverse process.

Let us call , , the -stable subordinator (see for details Bertoin [5], cap. III) and let us define its inverse process as its hitting time

 Eα(t)=inf{s>0:Vα(s)>t}. (2.15)

The processes and are characterized by their Laplace transforms. For the -stable subordinator we have

 Ee−ξVαt=e−tξα,α∈(0,1], (2.16)

and for its inverse process the time-Laplace transform reads

 ∫∞0e−ξt(Pr{Eα(t)∈ds}/ds)dt=ξα−1e−sξαα∈(0,1]. (2.17)
###### Theorem 2.2.

Let , , be the classical queue and let , , , be an inverse -stable subordinator (2.15) independent of . The fractional queue , , , can be represented as

 Nα(t)=N1(Eα(t)),t≥0,α∈(0,1), (2.18)

where the equality holds for the one-dimensional distribution.

###### Proof.

Let us consider the initial value problem

 {zDαtGα(z,t)=(1−z)[Gα(z,t)(μ−λz)−μpα0(t)],Gα(z,0)=zi, (2.19)

which is equivalent to (2.2). Applying the Laplace transform we obtain

 z[sα~Gα(z,s)−sα−1Gα(z,0)]=(1−z)[~Gα(z,s)(μ−λz)−μ~pα0(s)], (2.20)

Note that if (2.18) holds we can write

 ~Gα(z,s) =∫∞0e−st[∞∑k=0zk∫∞0Pr{Nα(y)=k}Pr{Eα(t)∈dy}]dt (2.21) =∫∞0e−st[∫∞0G(z,y)Pr{Eα(t)∈dy}]dt =∫∞0G(z,y)sα−1e−ysαdy,

and

 ~pα0(s) =∫∞0e−stpα0(t)dt (2.22) =∫∞0e−st[∫∞0p0(y)Pr{Eα(t)∈dy}]dt =∫∞0p0(y)sα−1e−ysαdy.

We now show that (2.21) and (2.22) satisfy (2.20). Observe that

 z[sα∫∞0G(z,y)e−ysαdyμ−zi]=(1−z)[(μ−λz)∫∞0G(z,y)e−ysαdy−μ∫∞0p0(y)e−ysαdy]. (2.23)

Consider the right hand side of (2.21). We can write

 (1−z)[(μ−λz)∫∞0G(z,y)e−ysαdy−μ∫∞0p0(y)e−ysαdy] (2.24) =∫∞0e−ysα(1−z)[(μ−λz)G(z,y)−μp0(y)]dy.

Considering that and satisfy

 z∂∂yG(z,y)=(1−z)[(μ−λz)G(z,y)−μp0(y)], (2.25)

we immediately obtain that

 (1−z)[(μ−λz)∫∞0G(z,y)e−ysαdy−μ∫∞0p0(y)e−ysαdy] (2.26) =z∫∞0e−ysα∂∂yG(z,y)dy =z[G(z,y)e−ysα∣∣y=∞y=0+sα∫∞0G(z,y)e−ysαdy] =z[sα∫∞0G(z,y)e−ysαdy−zi].

This concludes the proof. ∎

Using the Laplace transform (2.3) and the calculations carried out in Bailey [3] we can gain some insights on the mean value of the process.

###### Theorem 2.3.

We have that

 ENα(t)=i+(λ−μ)tαΓ(α+1)+μJαpα0(t), (2.27)

where

 Jαf(t)=1Γ(α)∫t0(y−t)α−1f(y)dy,t>0, (2.28)

is the Riemann–Liouville fractional integral [16].

###### Proof.

By means of the Laplace transform (2.3) of the probability generating function we can write

 ~ENα(s) =ddz~Gα(z,s)∣∣∣z=1 (2.29) =sα−1[−λ(z−a1)(z−a2)][(i+1)zi+ai+12(1−a2)−1]+λ(2−a1−a2)λ2(1−a1)2(1−a2)2 =sα−1[−i+1+ai+12(1−a2)−1λ(1−a1)(1−a2)+2−a1−a2λ(1−a1)2(1−a2)2] =sα−1[i+1+ai+12(1−a2)−1sα+2λ−(sα+λ+μ)λ2(1−a1)2(1−a2)2] =sα−1[i+1+ai+12(1−a−12)sα+2λ−(sα+λ+μ)s2α] =sα−1[1+i+ai+12(1−a2)−1sα+λ−μs2α−sαs2α] =sα−1[isα+λ−μs2α+ai+12(1−a2)−1sα] =is+λ−μsα+1+sα−1ai+12(1−a2)−1sα =is+λ−μsα+1+μ~pα0(s)sα.

Note the in the above calculation we have used the relation (2.10). Result (2.29) immediately yields

 ENα(t)=i+(λ−μ)tαΓ(α+1)+μJαpα0(t). (2.30)

###### Remark 2.1.

The validity of Theorem 2.2 can be checked with the aid of formula (2.27) as follows.

 ENα(t) =∫∞0EN1(w)Pr{Eα(t)∈dw} (2.31) =i+(λ−μ)∫∞0wPr{Eα(t)∈dw}+μ∫∞0∫t0p10(y)dyPr{Eα(t)∈dw}.

Therefore the time-Laplace transform, recalling that , can be written as

 ~ENα(s) =is+∫∞0wsα−1e−wsαdw+μ∫∞0∫w0p10(y)dysα−1e−wsαdw (2.32) =is+λ−μsα+1+μsα−1∫∞0p10(y)dy∫∞ye−wsαdw =is+λ−μsα+1+μsα−1∫∞0p10(y)dye−ysαsα =is+λ−μsα+1+μ~p10(sα)s.

Now, by noticing that (see formula (2.12)) we arrive at

 ~ENα(s)=is+λ−μsα+1+μ~pα0(s)sα, (2.33)

###### Remark 2.2.

A different form of formula (2.29) can be achieved by writing

 ~ENα(s) =is+λ−μsα+1+sα−1ai+12sα(1−a2) (2.34) =is+λ−μsα+1−sα−1λai+12(1−a1)s2α,

where we used the fact that . Furthermore, after considering

 a1=sα+λ+μλ−a2=μλa2, (2.35)

we arrive at

 ~ENα(s) =is+λ−μsα−1+sα−1ai2(μ−λa2)s2α (2.36) =is+λ−μsα−1+sα−1μai2−λai+12s2α =is+λ−μsα−1+μai2−λai+12sα+1.

Let us now address the problem of finding explicit results for the state probabilities of the proposed fractional queue model. We start by using the subordination relation stated in Theorem (2.2) with the classical solution of the M/M/1 queue in terms of modified Bessel functions of the first kind. In the non-fractional case () we have [3, Page 154]

 p1k(t)= (λμ)12(k−2)e−(λ+μ)tIi−k(2(λμ)1/2t) (2.37) +(λμ)12(k−i)∫t0e−(λ+μ)τ{λIi+k+2(2(λμ)1/2τ)−2(λμ)1/2Ii+k+1(2(λμ)1/2τ) +μIi+k(2(λμ)1/2τ)}dτ,

where is the modified Bessel function of the first kind.

The state probabilities , , , can thus be determined formally by subordination in the following way:

 pαk(t)=∫∞0p1k(y)Pr{Eα(t)∈dy}. (2.38)

Using the time-Laplace transform we have

 ~pαk(s)= ∫∞0p1k(y)sα−1e−ysαdy (2.39) = (λμ)12(k−2)∫∞0e−(λ+μ)yIi−k(2(λμ)1/2y)sα−1e−ysαdy +(λμ)12(k−i)∫∞0[∫y0e−(λ+μ)τ{λIi+k+2(2(λμ)1/2τ)−2(λμ)1/2Ii+k+1(2(λμ)1/2τ) +μIi+k(2(λμ)1/2τ)}dτ]sα−1e−ysαdy = (λμ)12(k−2)sα−1∫∞0e−y(sα+λ+μ)Ii−k(2(λμ)1/2y) +(λμ)12(k−i)sα−1∫∞0e−(λ+μ)τ{λIi+k+2(2(λμ)1/2τ)−2(λμ)1/2Ii+k+1(2(λμ)1/2τ) +μIi+k(2(λμ)1/2τ)}dτ∫∞τe−ysαdy,

Applying the well-known Laplace transform for we get

 ~pαk(s)= (λμ)12(k−2)sα−1[2(λμ)1/2]k−i[sα+λ+μ−√(sα+λ+μ)2−[2(λμ)1/2]2]i−k (2.40) ×[(sα+λ+μ)2−[2(λμ)1/2]2]−12 +(λμ)12(k−i)1s[λ∫∞0e−(sα+λ+μ)τIi+k+2(2(λμ)1/2τ)dτ −2(λμ)1/2∫∞0e−(sα+λ+μ)τIi+k+1(2(λμ)1/2τ) +μ∫∞0e−(sα+λ+μ)τIi+k(2(λμ)1/2τ)dτ] = (λμ)12(k−2)sα−1[sα+λ+μ−√(sα+λ+μ)2−4λμ2(λμ)1/2]i−k[(sα+λ+μ)2−4λμ]−1/2 +(λμ)12(k−i)λs[2(λμ)1/2]−(i+k+2)[sα+λ+μ−√(sα+λ+μ)2−4λμ]i+k+2 ×[(sα+λ+μ)2−4λμ]−1/2 −(λμ)12(k−i)2(λμ)1/2s[2(λμ)1/2]−(i+k+1)[sα+λ+μ−√(sα+λ+μ)2−4λμ]i+k+1 ×[(sα+λ+μ)2−4λμ]−1/2 +(λμ)12(k−i)μs[2(λμ)1/2]−(i+k)[sα+λ+μ−√(sα+λ+μ)2−4λμ]i+k ×[(sα+λ+μ)2−4λμ]−1/2 = sα−1(λ/μ)12(k−2)√(sα+λ+μ)2−4λμ[sα+λ+μ−√(sα+λ+μ)2−4λμ2(λμ)1/2]i−k +sα−1(λ/μ)12(k−1)λsα√(sα+λ+μ)2−4λμ[sα+λ+μ−√(sα+λ+μ)2−4λμ2(λμ)1/2]i+k+2 −sα−1(λ/μ)12(k−1)2(λμ)1/2sα√(sα+λ+μ)2−4λμ[sα+λ+μ−√(sα+λ+μ)2−4λμ2(λμ)1/2]i+k+1 +sα−1(λ/μ)12(k−1)μsα√(sα+λ+μ)2−4λμ[sα+λ+μ−√(sα+λ+μ)2−4λμ2(λμ)1/2]i+k.

Although the obtained Laplace transform has a clear structure it cannot be inverted in a simple manner. Note anyway that it should be related to the Laplace transform of some generalizations of Bessel functions.

In order to obtain more explicit results we must abandon the classical form of the state probabilities in terms of Bessel functions. We exploit instead a lesser known but certainly more appealing result due to Sharma [17, Chapter 2]. In particular we refer to equation (2.2.16) at page 17 which we recall here for the reader’s convenience. Here .

 p1k(t)= Missing or unrecognized delimiter for \right (2.41) +e−(λ+μ)t∞∑r=0(λt)k+r−i(μt)r(1r!(k+r−i)!−1(k+r)!(r−i)!).

By means of the above formula in the next theorem we derive an explicit expression for the state probabilities , , ,