Fractional Erlang Queues

TL;DR

This paper introduces a fractional generalization of Erlang queues using inverse stable subordinators, deriving equations, analyzing properties, and providing simulation algorithms for this new queue model.

Contribution

It presents the first fractional Erlang queue model, extending classical queues with fractional calculus and stable subordinators, and offers analytical and simulation tools.

Findings

Derived fractional Kolmogorov forward equations

Analyzed mean queue length and busy period distribution

Developed algorithms for simulating fractional queue paths

Abstract

We introduce a fractional generalization of the Erlang Queues $M/E_k/1$. Such process is obtained through a time-change via inverse stable subordinator of the classical queue process. We first exploit the (fractional) Kolmogorov forward equation for such process, then we use such equation to obtain an interpretation of this process in the queuing theory context. Then we also exploit the transient state probabilities and some features of this fractional queue model, such as the mean queue length, the distribution of the busy periods and some conditional distributions of the waiting times. Finally, we provide some algorithms to simulate their sample paths.

Equations

$$S=\{(n,s)\in N\times N:n\ge 1,1\le s\le k\}\cup \{(0,0)\}$$

$$p_{n,s}(t)$$

$$p_0(t)$$

$$\{\begin{array}{ll}\frac{d{p}_0}{dt}(t)=-\lambda p_0(t)+k\mu p_{1,1}(t)& \\ \frac{d{p}_{1,s}}{dt}(t)=-(\lambda +k\mu )p_{1,s}(t)+k\mu p_{1,s+1}(t)& 1\le s\le k-1\\ \frac{d{p}_{1,k}}{dt}(t)=-(\lambda +k\mu )p_{1,k}(t)+k\mu p_{2,1}(t)+\lambda p_0(t)& \\ \frac{d{p}_{n,s}}{dt}(t)=-(\lambda +k\mu )p_{n,s}(t)+k\mu p_{n,s+1}(t)+\lambda p_{n-1,s}(t)& n\ge 2,1\le s\le k-1\\ \frac{d{p}_{n,k}}{dt}(t)=-(\lambda +k\mu )p_{n,k}(t)+k\mu p_{n+1,1}(t)+\lambda p_{n-1,k}(t)& n\ge 2\\ p_0(0)=1& \\ p_{n,s}(0)=0& n\ge 1,1\le s\le k.\end{array}$$

$$I_n^k(t)={\left(\frac{t}{2}\right)}^n\sum _{r=0}^{+\infty }\frac{{\left(\frac{t}{2}\right)}^{r(k+1)}}{r!\Gamma (n+rk+1)},t\in C,n=0,1,2,\dots ,k=1,2,\dots$$

$$\begin{gathered} I_n^{k,s}(t)={\left(\frac{t}{2}\right)}^{n+k-s}\sum _{r=0}^{+\infty }\frac{{\left(\frac{t}{2}\right)}^{r(k+1)}}{(k(r+1)-s)!\Gamma (n+r+1)} \\ t\in C,n=0,1,2,\dots ,k=1,2,\dots ,s=1,2,\dots ,k. \end{gathered}$$

$$p_0(t)=\sum _{m=1}^{+\infty }m{\left(\frac{\lambda }{k\mu }\right)}^{-\frac{m}{k+1}}\frac{I_m^k\left(2{\left(\frac{\lambda }{k\mu }\right)}^{\frac{1}{k+1}}k\mu t\right)}{k\mu t}{e}^{-(\lambda +k\mu )t}$$

$$\begin{array}{rl}{p}_{n,s}(t)& ={\left(\frac{\lambda }{k\mu }\right)}^{\frac{k(n-1)+s}{k+1}}{I}_n^{k,s}(2(\lambda (k\mu {)}^k{)}^{\frac{1}{k+1}}t)e^{-(\lambda +k\mu )t}\\ & +k\mu {\left(\frac{\lambda }{k\mu }\right)}^{\frac{k(n-1)+s}{k+1}}{\int }_0^t{p}_0(z)I_n^{k,s}(2(\lambda (k\mu {)}^k{)}^{\frac{1}{k+1}}(t-z))e^{-(\lambda +k\mu )(t-z)}dz\\ & -k\mu {\left(\frac{\lambda }{k\mu }\right)}^{\frac{k(n-1)+s+1}{k+1}}{\int }_0^t{p}_0(z)I_n^{k,s+1}(2(\lambda (k\mu {)}^k{)}^{\frac{1}{k+1}}(t-z))e^{-(\lambda +k\mu )(t-z)}dz\end{array}$$

$$\begin{array}{rl}{p}_{n,k}(t)& ={\left(\frac{\lambda }{k\mu }\right)}^{\frac{kn}{k+1}}{I}_n^{k,k}(2(\lambda (k\mu {)}^k{)}^{\frac{1}{k+1}}t)e^{-(\lambda +k\mu )t}\\ & +k\mu {\left(\frac{\lambda }{k\mu }\right)}^{\frac{kn}{k+1}}{\int }_0^t{p}_0(z)I_n^{k,k}(2(\lambda (k\mu {)}^k{)}^{\frac{1}{k+1}}(t-z))e^{-(\lambda +k\mu )(t-z)}dz\\ & -k\mu {\left(\frac{\lambda }{k\mu }\right)}^{\frac{kn+1}{k+1}}{\int }_0^t{p}_0(z)I_{n+1}^{k,1}(2(\lambda (k\mu {)}^k{)}^{\frac{1}{k+1}}(t-z))e^{-(\lambda +k\mu )(t-z)}dz.\end{array}$$

$$p_0(t)=\sum _{m=1}^{+\infty }\sum _{r=0}^{+\infty }\frac{m{\lambda }^r(k\mu {)}^{m+rk-1}}{r!\Gamma (n+rk+1)}{t}^{m+r(k+1)-1}{e}^{-(\lambda +k\mu )t}$$

$$\begin{array}{rl}{p}_{n,s}(t)& =\sum _{j=0}^{+\infty }\frac{{\lambda }^{n+j}(k\mu {)}^{k(j+1)-s}}{(k(j+1)-s)!\Gamma (n+j+1)}{t}^{n+k-s+j(k+1)}{e}^{-(\lambda +k\mu )t}\\ & +\sum _{j=0}^{+\infty }\frac{{\lambda }^{n+j}(k\mu {)}^{k(j+1)-s+1}}{(k(j+1)-s)!\Gamma (n+j+1)}\\ & \times {\int }_0^t{p}_0(z)(t-z{)}^{n+k-s+j(k+1)}{e}^{-(\lambda +k\mu )(t-z)}dz\\ & -\sum _{j=0}^{+\infty }\frac{{\lambda }^{n+j}(k\mu {)}^{k(j+1)-s}}{(k(j+1)-s-1)!\Gamma (n+j+1)}\\ & \times {\int }_0^t{p}_0(z)(t-z{)}^{n+k-s-1+j(k+1)}{e}^{-(\lambda +k\mu )(t-z)}dz\end{array}$$

$$\begin{array}{rl}{p}_{n,k}(t)& =\sum _{j=0}^{+\infty }\frac{{\lambda }^{n+j}(k\mu {)}^{kj}}{(kj)!\Gamma (n+j+1)}{t}^{n+j(k+1)}{e}^{-(\lambda +k\mu )t}+\\ & +\sum _{j=0}^{+\infty }\frac{{\lambda }^{n+j}(k\mu {)}^{kj+1}}{(kj)!\Gamma (n+j+1)}\times {\int }_0^t{p}_0(z)(t-z{)}^{n+j(k+1)}{e}^{-(\lambda +k\mu )(t-z)}dz\\ & -\sum _{j=0}^{+\infty }\frac{{\lambda }^{n+j+1}(k\mu {)}^{k(j+1)}}{(k(j+1)-1)!\Gamma (n+j+2)}\times {\int }_0^t{p}_0(z)(t-z{)}^{n+k+j(k+1)}{e}^{-(\lambda +k\mu )(t-z)}dz.\end{array}$$

$$L(t)=\{\begin{array}{ll}k(N(t)-1)+S(t)& N(t)>0\\ 0& N(t)=0.\end{array}$$

$$m_k(n,s)=\{\begin{array}{ll}k(n-1)+s& (n,s)\in S^{\ast }\\ 0& (n,s)=(0,0)\end{array}$$

$$s_k(m)=\{\begin{array}{ll}\min \{s>0:s\in [m{]}_k\}& m>0\\ 0& m=0\end{array}$$

$$n_k(m)=\{\begin{array}{ll}\frac{m-s_k(m)}{k}+1& m>0\\ 0& m=0.\end{array}$$

$$P_n(t)=P(L(t)=n|L(0)=0).$$

$$\{\begin{array}{ll}\frac{d{P}_0}{dt}(t)=-\lambda P_0(t)+k\mu P_1(t)& \\ \frac{d{P}_n}{dt}(t)=-(\lambda +k\mu )P_n(t)+k\mu P_{n+1}(t)+\lambda {\sum }_{m=1}^n{c}_m{P}_{n-m}(t)& n\ge 1\\ P_n(0)={\delta }_{n,0}& n\ge 0\end{array}$$

$$P_m(t)=p_{n_k(m),s_k(m)}(t)$$

$$\{\begin{array}{l}\frac{dM}{dt}(t)=k(\lambda -\mu )+k\mu P_0(t)\\ M(0)=0\end{array}$$

$$M(t)=k(\lambda -\mu )t+k\mu {\int }_0^t{P}_0(y)dy.$$

$$B(t)=\sum _{r=0}^{+\infty }\frac{k{\lambda }^r(k\mu {)}^{k(r+1)}}{r!\Gamma (rk+k+1)}{\int }_0^t{z}^{k+r(k+1)-1}{e}^{-(\lambda +k\mu )z}dz.$$

$$w(\xi ;t)=\sum _{n=0}^{+\infty }{P}_n(t)k\mu \frac{(k\mu t{)}^{n-1}}{(n-1)!}{e}^{-k\mu t}.$$

$$E[e^{-v{\sigma }_{\nu }(t)}]=e^{-t{v}^{\nu }},v>0,\nu \in (0,1).$$

$$L_{\nu }(t)=\inf \{s\ge 0:{\sigma }_{\nu }(s)>t\}t\ge 0$$

$$Q^{\nu }(t):=(N^{\nu }(t),S^{\nu }(t))=Q(L_{\nu }(t)).$$

$$p_{n,s}^{\nu }(t)$$

$$p_0^{\nu }(t)$$

$$I_t^{\nu }x=\frac{1}{\Gamma (\nu )}{\int }_0^t(t-\tau {)}^{\nu -1}x(\tau )d\tau$$

$$D_t^{\nu }x=\frac{1}{\Gamma (1-\nu )}{\int }_0^t\frac{dx}{dt}(\tau )\frac{d\tau }{(t-\tau {)}^{\nu }}.$$

$$L[D_t^{\nu }x](z)=z^{\nu }\overline{x}(z)-z^{\nu -1}x(0).$$

Full text

N. Leonenko was supported in particular by Australian Research Council’s Discovery Projects funding scheme (project DP160101366), and by project MTM2015-71839-P of MINECO, Spain (co-funded with FEDER funds).

1. Introduction

The Erlang Queue is one of the most popular model for a queue process. It is characterized by arrivals determined by a Poisson process and independent Erlang service times.

We introduce a non-Markovian generalization of the classical Erlang $M/E_{k}/1$ queue process based on the fractional Poisson process, FPP for short ([1]). There are essentially three approaches to the concept of FPP. The “renewal” approach consists of considering the characterization of the Poisson process as a sum of independent non-negative random variables and, instea of the assumption that these random variables have an exponential distribution, we assume that they have the Mittag-Leffler distribution ([22, 26, 27]). In [4, 5], the renewal approach to the FPP is developed and it is proven that its $n$-dimensional distributions coincide with the solution to “fractionalized” differential-difference equations with fractional derivatives in time. Finally, using “inverse subordination”, a FPP can be constructed as a time-changed classical Poisson process ([28]).

The first fractional generalization of the classical $M/M/1$ queue process was proposed in [7]. The transient behaviour of the fractional $M/M/1$ queues with catastrophes (which in the classical case are studied for instance in [10] and further generalized in [8, 9, 14]) is investigated in [2].

In this paper, we present a study of the transient behaviour of the fractional Erlang queue $M/E_{k}/1$.

In Section 2, we recall the basic definitions and properties for a classical Erlang queue $M/E_{k}/1$. In particular we recall the formulas for the transient state probabilities, the mean queue length and the distribution of the busy period.

In Section 3 we give the definition of Fractional Erlang Queue and we determine the fractional forward Kolmogorov equation for the transient state probabilities of such queue.

In Section 4 the Fractional Erlang Queue process is investigated in terms of inter-arrival times and service times. Such investigation leads to an interpretation of this process in a queueing theory context and the characterization of the distributions of the inter-arrival, inter-phase and service times.

In Section 5 formulas for transient state probabilities of such process and their Laplace transform are finally obtained. In Section 6 some characteristics of the process are studied. In particular in Subsection 6.1, a fractional differential equation for the mean queue length is provided. Moreover, we find also a formula for the mean queue length and its Laplace transform and, by using such formula, we provide an expression for the fractional integral of the probability that the queue is empty. In Subsection 6.2 we provide also a formula for the distribution of the busy period. In Subsection 6.3 we focus on the waiting times, determining a formula for a particular conditional waiting time.

Finally, in Section 7, we provide an algorithm for sample paths simulation of such process, which is a modified version of the well-known Gillespie algorithm (see [13]).

2. Erlang Queues

Let us recall the definition of Erlang queue. A $M/E_{k}/1$ queue is a queue with Poisson input (and so exponential independent interarrival times) with parameter $\lambda$ and an Erlang service system with shape parameter $k$ and rate $\mu$. Such service system can be seen as a system composed of $k$ phases each with exponential service time with parameter $k\mu$.

Let us denote with $N(t)$ the number of the customers in the system at the time $t\geq 0$. Moreover, if the queue is not empty, let us denote with $S(t)$ the phase of the customer that is currently being served at the time $t\geq 0$. If $N(t)=0$ let us denote also $S(t)=0$. Finally denote $Q(t)=(N(t),S(t))$ for $t\geq 0$, that is a sort of state-phase process with state space:

[TABLE]

and let us denote $\operatorname{\mathcal{S}}^{*}=\{(n,s)\in\operatorname{\mathbb{N}}\times\operatorname{\mathbb{N}}:\ n\geq 1,\ 1\leq s\leq k\}$. Let us define the transient state-phase probability functions as

[TABLE]

It is well known (see for instance [34]) that the function $p_{n,s}(t)$ and $p_{0}(t)$ satisfy the following difference-differential equations:

[TABLE]

In [16] a solution of this system is obtained by means of generalized modified Bessel functions defined in [25] as

[TABLE]

and in [15] as

[TABLE]

In particular we have:

[TABLE]

while for $n\geq 1$ and $1\leq s\leq k-1$:

[TABLE]

and finally for $n\geq 1$ and $s=k$:

[TABLE]

It will be useful in the following to write the Bessel functions explicitly as power series. In this way we have

[TABLE]

while for $n\geq 1$ and $1\leq s\leq k-1$ we have

[TABLE]

and finally for $n\geq 1$ and $s=k$ we have

[TABLE]

We can also give an equivalent representation of the queue by using the queue length process, that is to say the process $\operatorname{\mathcal{L}}(t)$ given by:

[TABLE]

The queue length process gives us the length of the queue by means of phases. In particular, since the map $m_{k}:\operatorname{\mathcal{S}}\to\operatorname{\mathbb{N}}_{0}$ given by:

[TABLE]

is bijective, each state of $\operatorname{\mathcal{L}}(t)$ represents only one state of $Q(t)$ and vice-versa, so that we can equivalently describe the queue with $Q(t)$ or $\operatorname{\mathcal{L}}(t)$: indeed we have $m_{k}(Q(t))=\operatorname{\mathcal{L}}(t)$. Let us also recall that the inverse map of $m_{k}(n,s)$ is given by $(n_{k}(m),s_{k}(m))$ where:

[TABLE]

where $[m]_{k}$ is the residual class of $m$ modulo $k$ and:

[TABLE]

Let us define the state probabilities of $\operatorname{\mathcal{L}}(t)$ as

[TABLE]

In [24] and [25] the forward equations for $\operatorname{\mathcal{P}}_{n}$ are given as:

[TABLE]

where $c_{m}=\delta_{m,k}$ for $M/E_{k}/1$ queues. However, the bijectivity of $m_{k}(n,s)$ gives us the equality

[TABLE]

so that the difference-differential system (11) is just a compact rewriting of (1) where the states of $Q(t)$ are disposed in lexicographic order.

In [24] the mean queue length, defined as $M(t)=\operatorname{\mathbb{E}}[\operatorname{\mathcal{L}}(t)|\operatorname{\mathcal{L}}(0)=0]$, is show to solve the following Cauchy problem:

[TABLE]

which leads, in [25], to the formula:

[TABLE]

Moreover, in [25], the distribution function of the busy period has been obtained as

[TABLE]

Finally, let us also recall that the probability density function of the waiting time of a customer that enters the system at time $t>0$ has been given in [11] as:

[TABLE]

The author obtained such formula by using explicitly the Markov property of the queue.

3. Fractional Erlang Queues

In [7] the authors introduce a fractional version of the $M/M/1$ queue. In particular, such queue has been shown useful to model some financial data for which classical birth-death processes were not enough. Following this idea, now we give the definition of a fractional $M/E_{k}/1$ queue. Let us consider a classical state-phase process $Q(t)=(N(t),S(t))$ for the $M/E_{k}/1$ queue and define $\{\sigma_{\nu}(t)\}_{t\geq 0}$ a $\nu$-stable subordinator independent from $Q(t)$, such that (see, for instace, [29])

[TABLE]

Define

[TABLE]

the inverse of the $\nu$-stable subordinator and pose $N^{\nu}(t):=N(L_{\nu}(t))$, $S^{\nu}(t):=S(L_{\nu}(t))$ and

[TABLE]

We will say that the process $Q^{\nu}(t)$ is the state-phase process of the fractional Erlang queue $M/E_{k}/1$. By definition, let us also pose $Q^{1}(t)=Q(t)$.

Before giving an interpretation in the sense of the queueing theory, let us show what are the forward equations for the state probabilities of $Q^{\nu}(t)$. Let us define:

[TABLE]

To obtain such forward equations, we will need the fractional integral (see [23]) of order $\nu\in(0,1)$ defined, for a function $x:[0,t_{1}]\subseteq\operatorname{\mathbb{R}}\to\operatorname{\mathbb{R}}$, as

[TABLE]

for any $t\in[0,t_{1}]$ and the Caputo fractional derivative of order $\nu\in(0,1)$

[TABLE]

The classes of functions for which such operators are well-defined are discussed in [29].

Moreover, we will use a Laplace transform approach, as described in [19, 33]. Let us also remark that if $\overline{x}$ is the Laplace transform of the function $x$ and $\operatorname{\mathcal{L}}$ is the Laplace transform operator, then:

[TABLE]

For the state probabilities we can show the following theorem.

Theorem 3.1**.**

The state probabilities $(p^{\nu}_{n,s})_{(n,s)\in\operatorname{\mathcal{S}}}$ are solution of the following fractional difference-differential Cauchy problem:

[TABLE]

where $D_{t}^{\nu}$ is the Caputo fractional derivative.

Proof.

Let us define the probability generating function

[TABLE]

By multiplying the second equation of (18) by $z^{s+1}$, the third equation of (18) by $z^{k+1}$, the fourth equation of (18) by $z^{k(n-1)+s+1}$, the fifth equation of (18) by $z^{kn+1}$ and the summing all these equations with respect to $n$ and $s$ we have that the state probabilities $p^{\nu}_{n,s}$ solve this Cauchy problem if and only if $G^{\nu}(z,t)$ solves the following Cauchy problem:

[TABLE]

Denoting with $\widetilde{G}^{\nu}(z,v)$ and $\pi_{0}^{\nu}(v)$ the Laplace transform of $G^{\nu}(z,t)$ and $p_{0}^{\nu}(t)$, we know that $G^{\nu}(z,t)$ solves the Cauchy problem (19) if and only if its Laplace transform solves:

[TABLE]

where we used the formula for the Laplace transform of the Caputo derivative given in formula (17). Now let us observe that:

[TABLE]

and, analogously:

[TABLE]

Thus for the probability generating function we know that

[TABLE]

Let us recall from [30] the following Laplace transform

[TABLE]

Thus by using Eq. (24) in Eqs. (22) and (23) we obtain:

[TABLE]

Using these two formulas we have that $\widetilde{G}^{\nu}(z,v)$ solves Eq. (20) if and only if

[TABLE]

From [16] we know that $G^{1}(z,t)$ is solution of the following Cauchy problem:

[TABLE]

Using such property in (26) we obtain

[TABLE]

which, by using the integration by parts formula, gives us an identity. ∎

Let us introduce, as we did for the classical Erlang queue, the queue length process:

[TABLE]

and define $\operatorname{\mathcal{P}}_{m}^{\nu}(t):=\operatorname{\mathbb{P}}(\operatorname{\mathcal{L}}^{\nu}(t)=m|\operatorname{\mathcal{L}}^{\nu}(0)=0)$. By the equality $\operatorname{\mathcal{P}}_{m}^{\nu}(t)=p^{\nu}_{n_{k}(m),s_{k}(m)}(t)$ and Theorem 3.1, we know that $(\operatorname{\mathcal{P}}_{m}^{\nu}(t))_{m\in\operatorname{\mathbb{N}}}$ is solution of the fractional Cauchy problem

[TABLE]

where $c_{j}=\delta_{j,k}$, which is equivalent to the fractional Cauchy problem (18). Since the solutions $\operatorname{\mathcal{P}}^{\nu}_{m}(t)$ are defined as probabilities, by using Corollary $2$ of [2] and Schur’s test (see [17]) one can prove the following:

Corollary 3.2**.**

The fractional Cauchy problems (18) and (28) admit unique global solutions $\mathbf{p}^{\nu}(t)=(p_{n,s}^{\nu}(t))_{(n,s)\in\operatorname{\mathcal{S}}}$ and $\mathbf{P}^{\nu}(t)=(\operatorname{\mathcal{P}}_{m}^{\nu}(t))_{m\in\operatorname{\mathbb{N}}}$, where $\mathbf{p}^{\nu},\mathbf{P}^{\nu}:[0,+\infty)\to l^{2}$.

4. Interpretation of the fractional $M/E_{k}/1$ queue

By using the result given in Theorem 3.1 we can show what are the main features of a fractional $M/E_{k}/1$ queue in terms of interarrival and service times. Let us define:

[TABLE]

Note that

[TABLE]

Let us work with the functions $h_{n}(t)$. By summing the equations of (18) with respect to $s$ we obtain the following fractional Cauchy problem

[TABLE]

that, for $\nu=1$, gives us

[TABLE]

Before proceeding, let us give some notation. Let $T$ be a random variable with distribution function given by

[TABLE]

where $E_{\nu}$ is the Mittag-Leffler function (see [20]) defined as

[TABLE]

we call $T$ a Mittag-Leffler random variable with fractional index $\nu$ and parameter $\lambda$ and we denote it by $T\sim ML_{\nu}(\lambda)$. In particular let us denote $\Psi(t)=E_{\nu}(-\lambda t^{\nu})$, $f_{T}$ the probability density function of $T$ and denoting with $\overline{f}_{T}$ its Laplace transform and recalling the Laplace transform of the Mittag-Leffler function that can be obtained from (35) by posing $\delta=\gamma=1$, we have

[TABLE]

Now we can show the following Theorem.

Theorem 4.1**.**

The interarrival times of a $\nu$-fractional $M/E_{k}/1$ queue are i.i.d. Mittag-Leffler random variables with fractional index $\nu$ and parameter $\lambda$.

Proof.

The independence of the interarrival times of a fractional $M/E_{k}/1$ queue follows easily from the definition. Let us show that these interarrival times are Mittag-Leffler distributed. To do this, let us introduce the process $\widetilde{N}^{1}(t)$ as a modification of $N^{1}(t)$ such that the state probabilities

[TABLE]

are solution of the following Cauchy problem

[TABLE]

that is to say that $\widetilde{N}^{1}(t)$ is a birth process ($\mu=0$) that has the same birth rate $\lambda$ as $N^{1}(t)$, starts at $m$ ad admits $m+1$ as absorbent state. In particular the Cauchy problem can be reduced to a couple of linear ODEs:

[TABLE]

Moreover, since they are probabilities, we have also that $\widetilde{h}_{m}^{1}(t)+\widetilde{h}_{m+1}^{1}(t)=1$ for any $t\geq 0$. Finally, denoting with $T^{1}$ the arrival time of the first new customer of the queue $Q^{1}$ starting from $N^{1}(0)=m$, we have that $\widetilde{h}_{m+1}^{1}(t)=F_{T^{1}}(t)$.

Now let us define $\widetilde{N}^{\nu}(t)=\widetilde{N}^{1}(L_{\nu}(t))$ and

[TABLE]

By definition, $\widetilde{N}^{\nu}(t)$ is a process that starts from $m$ and admits $m+1$ as absorbent state. Moreover, denoting with $T^{\nu}$ the arrival time of the first new customer in the queue $Q^{\nu}$ starting from $N^{\nu}(0)=m$, we have that $\widetilde{h}_{m+1}^{\nu}(t)=F_{T^{\nu}}(t)$. Indeed, since $N^{1}(t)$ is a Markov process, $N^{\nu}(t)$ is a semi-Markov process and, denoting with $T_{n}$ the $n$-th jump time, the discrete time process $N^{\nu}(T_{n})$ is a time-homogeneous Markov chain (see, for instance, [12]). For this reason, $T^{\nu}$ can be seen as the interarrival time of the queue $Q^{\nu}$.

Now we have to show that

[TABLE]

The initial conditions are obviously verified. Now let us work with the first equation. We know that this equation is verified if and only if, denoting with $\overline{h}_{m}^{\nu}$ the Laplace transform of $\widetilde{h}_{m}^{\nu}$:

[TABLE]

By using the definition of $\widetilde{h}_{m}^{\nu}(t)$ we have that:

[TABLE]

and thus, taking the Laplace transform

[TABLE]

Using this formula, we know that Eq. (32) is verified if and only if

[TABLE]

Recalling that $\widetilde{h}^{1}_{m}(t)$ is solution of the following Cauchy problem

[TABLE]

we have that Eq. (33) is verified if and only if

[TABLE]

that, by using the integration by parts formula, gives us an identity. In an analogous way, one can show that the second equation of (31) is verified.

Finally, we only have to solve the fractional Cauchy problem (31). To do this, let us observe that the first equation is independent from the second one, so we can solve first the Cauchy problem:

[TABLE]

whose solution is $\widetilde{h}_{m}^{\nu}(t)=E_{\nu}(-\lambda t^{\nu})$ (see [20]). Finally:

[TABLE]

∎

Now let us work with the functions $q_{s}^{\nu}(t)$, in order to obtain some information on the service times. Summing the equations in (18) with respect to $n$ we have:

[TABLE]

that, for $\nu=1$, gives us

[TABLE]

We can use the functions $q_{s}^{\nu}$ exactly in the same way we did with $h_{n}^{\nu}$ to show the following Theorem.

Theorem 4.2**.**

The inter-phase times of a $\nu$-fractional $M/E_{k}/1$ queue are i.i.d. Mittag-Leffler random variables with fractional index $\nu$ and parameter $k\mu$.

Let us give some other notation. Let $T$ be a random variable with distribution function given by

[TABLE]

we call $T$ a generalized Erlang random variable (as introduced in [26]) with fractional index $\nu$, shape parameter $k$ and rate $\lambda$ and we denote it by $T\sim GE_{\nu}(k,\lambda)$. In particular $T\sim GE_{\nu}(k,\lambda)$ is sum of $k$ i.i.d. Mittag-Leffler random variable of fractional index $\nu$ and parameter $\lambda$. For this reason, denoting with $f_{T}$ the probability density function of $T$ and $\overline{f}_{T}$ its Laplace transform, we have, from (30),

[TABLE]

From Theorem 4.2 we easily obtain the following Corollary

Corollary 4.3**.**

The service times of a $\nu$-fractional $M/E_{k}/1$ queue are i.i.d. generalized Erlang random variables with fractional index $\nu$, shape parameter $k$ and rate $k\mu$.

Proof.

Let us simply observe that a service time is sum of $k$ inter-phase times, which are i.i.d. Mittag-Leffler random variables of fractional index $\nu$ and parameter $k\mu$. ∎

From Theorems 4.1 and 4.2 and Corollary 4.3 we can conclude that a fractional $M/E_{k}/1$ queue is a single-channel queue with Mittag-Leffler interarrival times and generalized Erlang service times in which the service is given through $k$ phases, each one with Mittag-Leffler distributed service time. In particular for $k=1$ we obtain the fractional $M/M/1$ queue introduced in [7].

We can also obtain the distribution of a general sojourn time of $\operatorname{\mathcal{L}}^{\nu}(t)$ working with $\operatorname{\mathcal{P}}_{n}^{\nu}(t)$ as we already did with $q_{s}^{\nu}(t)$ and $h_{n}^{\nu}(t)$.

Theorem 4.4**.**

The sojourn time of the queue length process of a $\nu$-fractional $M/E_{k}/1$ queue in a state $\overline{n}>0$ is a Mittag-Leffler random variable with fractional index $\nu$ and parameter $\lambda+k\mu$.

Remark 4.5*.*

The sojourn time of $\operatorname{\mathcal{L}}^{\nu}(t)$ in the state [math] coincides with an interarrival time, so it is a Mittag-Leffler random variable of fractional index $\nu$ and parameter $\lambda$.

From this information on sojourn time, we can also obtain a new information on the mutual dependence of interarrival times and inter-phase times.

Corollary 4.6**.**

Interarrival times and inter-phase times are not independent for $\nu\in(0,1)$.

Proof.

Let us first remark that if two random variables $T$ and $S$ are independent, then

[TABLE]

Let $T_{S}$ be a sojourn time, $T_{A}$ an interarrival time and $T_{P}$ an inter-phase time. Thus we have $T_{S}=\min\{T_{A},T_{P}\}$. However

[TABLE]

thus we have that $\operatorname{\mathbb{P}}(T_{S}\geq t)=\operatorname{\mathbb{P}}(T_{P}\geq t)\operatorname{\mathbb{P}}(T_{A}\geq t)$ if and only if

[TABLE]

that is true if and only if $\lambda=k\mu=0$ or $\nu=1$ (see [32]). ∎

5. Transient state probabilities

Our aim is now to obtain the transient state probability functions $p_{0}^{\nu}(t)$ and $p_{n,s}^{\nu}(t)$ in a closed form, with the aid of the three-parameter Mittag-Leffler function $E_{\nu,\mu}^{\delta}$ (see [20]) defined as

[TABLE]

and the Laplace transform formula (see [18])

[TABLE]

To do this, let us first work with $p_{0}^{\nu}(t)$.

Theorem 5.1**.**

Let $p_{0}^{\nu}(t)$ be the [math] state probability of a fractional $M/E_{k}/1$ queue $Q^{\nu}(t)$ and $\pi_{0}^{\nu}(v)$ its Laplace transform. Then we have

[TABLE]

and

[TABLE]

where

[TABLE]

Proof.

Let us recall that we already have a closed form for $p_{0}^{1}(t)$, given by Eq. (5). Moreover, using this closed form in Eq. (25) we have:

[TABLE]

that, integrating, becomes

[TABLE]

Now let us pose $C^{0}_{m,r}$, $\delta^{0}_{m,r}$ and $\gamma^{0}_{m,r}$ as in (38) and observe that

[TABLE]

Finally, by using Eq. (35) we obtain:

[TABLE]

∎

Remark 5.2*.*

For $k=1$ we have that:

[TABLE]

and then we have

[TABLE]

Now, let us change variables, posing $\overline{m}=m+r$, so we have

[TABLE]

This is the second formula obtained in Remark $4$ in [2] for the fractional $M/M/1$ queue.

Remark 5.3*.*

For $\nu=1$ we obtain

[TABLE]

and then, with some calculations

[TABLE]

Let us remark that

[TABLE]

so that

[TABLE]

that is Eq. (5).

To work with $p_{n,s}^{\nu}(t)$, it could be useful to show the following Lemma.

Lemma 5.4**.**

For any $n\in\operatorname{\mathbb{N}}$ we have

[TABLE]

where $C^{0}_{m,r}$ is defined in (38).

Proof.

Let us first notice that:

[TABLE]

and then we can use Fubini’s theorem to obtain

[TABLE]

Now let us pose $w=y-z$ to obtain

[TABLE]

that, since $\chi_{[0,z+w]}(z)=\chi_{[0,+\infty)}(w)$, gives us

[TABLE]

By integrating, we have

[TABLE]

Moreover, by (25) and (37) we have

[TABLE]

so finally we obtain

[TABLE]

∎

Now we are ready to show the following Theorem.

Theorem 5.5**.**

For any $(n,s)\in\operatorname{\mathcal{S}}^{*}$ let $p_{n,s}^{\nu}(t)=\operatorname{\mathbb{P}}(Q^{\nu}(t)=(n,s)|Q^{\nu}(0)=(0,0))$ and $\pi_{n,s}^{\nu}(v)$ be its Laplace transform. Then

[TABLE]

and

[TABLE]

where

[TABLE]

Proof.

From Eq. (21) we easily obtain

[TABLE]

Let us first consider $1\leq s\leq k-1$; by using Eq. (6), we have

[TABLE]

Integrating, using Eq. (39) and posing $A^{n,s}_{j}$, $a^{n,s}_{j}$, $\alpha^{n,s}_{j}$, $B^{n,s}_{j,m,r}$, $b^{n,s}_{j,m,r}$, $\beta^{n,s}_{j,m,r}$, $C^{n,s}_{j,m,r}$, $c^{n,s}_{j,m,r}$ and $\gamma^{n,s}_{j,m,r}$ as in equation (40) for $s\not=k$ we obtain

[TABLE]

Finally, by using Eq. (35) we obtain

[TABLE]

Now, for $s=k$, from Eq. (41) and (7) we obtain

[TABLE]

Integrating, using Eq. (39) and posing $A^{n,s}_{j}$, $a^{n,s}_{j}$, $\alpha^{n,s}_{j}$, $B^{n,s}_{j,m,r}$, $b^{n,s}_{j,m,r}$, $\beta^{n,s}_{j,m,r}$, $C^{n,s}_{j,m,r}$, $c^{n,s}_{j,m,r}$ and $\gamma^{n,s}_{j,m,r}$ as in equation (40) for $s=k$ we obtain

[TABLE]

Finally, by using Eq. (35) we obtain

[TABLE]

∎

6. Some features of the fractional $M/E_{k}/1$ queue

In this section we want to determine some features of the fractional $M/E_{k}/1$ queue. In particular we focus on the mean queue length (measured in phases), the distribution of the busy period and the distribution of the waiting times.

6.1. Mean queue length

To determine the mean queue length, let us consider the queue length process $\operatorname{\mathcal{L}}^{\nu}(t)$ and define:

[TABLE]

Let us first show the following Theorem

Theorem 6.1**.**

The function $M^{\nu}(t)$ is solution of the following fractional Cauchy problem:

[TABLE]

Proof.

The initial condition follows obviously from the definition of $M^{\nu}(t)$.

To obtain the equation, recall that, by definition:

[TABLE]

Thus, let us consider Eq. (28), multiply the second equation by $n$ and then sum over $n$ to obtain

[TABLE]

Let us observe that

[TABLE]

and

[TABLE]

Now let us pose $r=n-m$ to obtain

[TABLE]

Recalling that $c_{m}=\delta_{m,k}$, we have

[TABLE]

Using Eqs. (44) and (45) in (43) we obtain the first equation of (42). ∎

By integrating (42), we obtain the following representation of the mean queue length:

[TABLE]

For $k=1$ we obtain the mean queue length of a fractional $M/M/1$ queue, as stated in Theorem $2.3$ of [7]. For $\nu=1$ we obtain Eq. (12).

Starting from Eq. (12) we can also obtain a closed form for $M^{\nu}(t)$.

Theorem 6.2**.**

Let $M^{\nu}(t)=\operatorname{\mathbb{E}}[\operatorname{\mathcal{L}}^{\nu}(t)|\operatorname{\mathcal{L}}^{\nu}(0)=0]$ and denote with $\overline{M}^{\nu}(v)$ its Laplace transform. Then

[TABLE]

and

[TABLE]

where

[TABLE]

and $C^{0}_{m,r}$ and $\delta^{0}_{m,r}$ are defined in (38).

Proof.

From Eq. (21) and the equality $\operatorname{\mathcal{P}}_{m}^{\nu}=p_{n_{k}(m),s_{k}(m)}^{\nu}$ we easily obtain that

[TABLE]

Using (12) we obtain

[TABLE]

Taking the Laplace transform we obtain

[TABLE]

For the first integral we easily obtain

[TABLE]

For the second one we have

[TABLE]

By using Eq. (5) we obtain

[TABLE]

Using Eqs. (52), (51) and (50) in (49) we obtain

[TABLE]

Posing then $\gamma^{M}_{m,r}$ as in (48) we obtain

[TABLE]

Finally, taking the inverse Laplace transform and recalling eq. (35) we obtain

[TABLE]

∎

Comparing Eq. (46) with (47) we also obtain some information on the fractional integral of $\operatorname{\mathcal{P}}_{0}^{\nu}(t)$. In particular

[TABLE]

6.2. Distribution of the busy period

We can obtain the distribution of the busy period for a fractional $M/E_{k}/1$ queue by following the lines of [24].

Theorem 6.3**.**

Let $B^{\nu}$ be the duration of the busy period of a fractional $M/E_{k}/1$ queue, $B^{\nu}(t)$ be its distribution function and $\overline{B}^{\nu}(v)$ be its Laplace transform. Then we have

[TABLE]

and

[TABLE]

where

[TABLE]

with $C^{0}_{k,r}$ and $\delta^{0}_{k,r}$ defined in Eqs. (38).

Proof.

Let us consider the process $\widetilde{\operatorname{\mathcal{L}}}^{1}(t)$ whose state probabilities

[TABLE]

satisfy the following Cauchy problem:

[TABLE]

The process $\operatorname{\widetilde{\mathcal{L}}}^{1}(t)$behaves as $\operatorname{\mathcal{L}}^{1}(t)$ except that it starts from $k$ (that is to say when the first customer entered the queue) instead that from [math] and admits [math] as absorbent state. In particular, by construction, $\operatorname{\widetilde{\mathcal{P}}}_{0}^{1}(t)=B^{1}(t)$.

Now let us consider $\operatorname{\widetilde{\mathcal{L}}}^{\nu}(t):=\operatorname{\widetilde{\mathcal{L}}}^{1}(L_{\nu}(t))$. As done with the interarrival and the inter-phase times, if $T_{n}$ is the $n$-th jump time of the process $\operatorname{\mathcal{L}}^{\nu}(t)$, then $\operatorname{\mathcal{L}}^{\nu}(T_{n})$ is a time-homogeneous Markov chain, then, by construction, $\operatorname{\widetilde{\mathcal{P}}}_{0}^{\nu}(t)=B^{\nu}(t)$.

From the relation between $\operatorname{\widetilde{\mathcal{L}}}^{1}$ and $\operatorname{\widetilde{\mathcal{L}}}^{\nu}$, we have that

[TABLE]

Thus, by using Eq. (13) we have

[TABLE]

Taking the Laplace transform we have

[TABLE]

Now let us consider the integral and observe that:

[TABLE]

Using Eq. (55) in Eq. (54) we obtain

[TABLE]

Thus, posing $C^{B}_{r}$ and $\gamma^{B}_{r}$ as in (53) we obtain

[TABLE]

Finally, taking the inverse Laplace transform and using Eq. (35), we obtain

[TABLE]

∎

6.3. Distribution of some conditioned waiting times

For classical single-channel queues, a method to obtain the waiting time distribution has been presented in [11]. However, this technique makes an explicit use of the Markov property of a classical single-channel queue. Our fractional queues are semi-Markov processes, then this technique cannot be used in our context. However, to give some information on the waiting times we can introduce some form of conditioning.

Before working with the waiting times, let us introduce a new distribution function that will be useful in the following. Consider a random variable $T\sim ML_{\nu}(\lambda)$ and fix $t_{0}>0$. Observe that

[TABLE]

The explicit dependence of this probability with respect to $t_{0}$ is due to the lack of semigroup property that is typical of Mittag-Leffler functions (see, for instance, [32]). This observation leads us to the definition of a new distribution.

Consider a random variable $T$ such that its distribution function is given by

[TABLE]

we call such random variable residual Mittag-Leffler random variable starting from $t_{0}$ with fractional index $\nu$ and parameter $\lambda$ and we denote it by $T\sim RML_{\nu}(t_{0},\lambda)$. Denote

[TABLE]

and observe that, by the definition of Mittag-Leffler function given in (29), we have

[TABLE]

Denote with $\overline{\Phi}(v)$ its Laplace transform and observe that

[TABLE]

where $\Gamma(z,s)$ is the analytic extension of the upper incomplete gamma function on $\operatorname{\mathbb{C}}^{2}$ (see, for instance, [21, 31]). Denote with $f_{T}(t)$ the probability density function of $T$ and observe that

[TABLE]

Denoting with $\overline{f}_{T}(v)$ its Laplace transform, we obtain

[TABLE]

Finally let us also observe that if $T\sim RML_{\nu}(t_{0},\lambda)$ for $\nu=1$, then $T\sim Exp(\lambda)$. Now we are ready to determine such waiting times.

Let $W^{\nu}_{t}$ be the waiting time of a customer that enters the queue at time $t>0$. Let us define the random set

[TABLE]

that is the set of the time instants that precede $t$ in which a phase has been completed. Let us define also the random variable

[TABLE]

which is the number of phase terminated before $t$. Finally, let us define the random variable

[TABLE]

that is the last instant before $t$ in which a phase has been completed. Let us define:

[TABLE]

Since we have conditioned this probability with $\operatorname{\mathcal{L}}^{\nu}(t)=n$, we know that a customer has to wait $n$ independent phases to be terminated. Thus we have

[TABLE]

where $T_{j}$ are $n$ independent random variables. For $1\leq j\leq n-1$, these are just Mittag-Leffler random variables with $T_{j}\sim ML_{\nu}(k\mu)$, thus, denoting $\widetilde{W}=\sum_{j=1}^{n-1}T_{j}$, we know that $\widetilde{W}\sim GE_{\nu}(n-1,k\mu)$. For $T_{n}$, we know that the last time a phase has been completed is $t_{0}$, since we have conditioned the probability with $\Theta_{t}^{\nu}=t_{0}$. So we know that a customer is in the current phase for $t-t_{0}$. For such reason, we know that $T_{n}\sim RML_{\nu}(t-t_{0},k\mu)$. Thus we have that

[TABLE]

In particular we can obtain the Laplace transform $\overline{w}^{\nu}(v;t,t_{0},n)$ of $w^{\nu}(\xi;t,t_{0},n)$ as

[TABLE]

For $\nu=1$, we know that $T_{n}\sim Exp(k\mu)$ and $\widetilde{W}\sim E_{n-1}(k\mu)$, thus $T_{n}+\widetilde{W}\sim E_{n}(k\mu)$ and we know that

[TABLE]

since the density is independent from $t_{0}$. Thus we obtain

[TABLE]

that is Eq. (14), that is the one given in [11].

7. Sample paths simulation

In this section we want to give some algorithm to simulate such queue process. To do this, we need to show how to simulate a Mittag-Leffler random variable. Let us show two preliminary simple lemmas.

Lemma 7.1**.**

Let $T\sim Exp(\lambda)$ and $\sigma_{\nu}(t)$ be a $\nu$-stable subordinator independent of $T$. Then $\sigma_{\nu}(T)\sim ML_{\nu}(\lambda)$

Proof.

Let us just observe that

[TABLE]

where the last equality is given by the Laplace transform of the density $\operatorname{\mathbb{P}}(L_{\nu}(t)\in dy)$ of the inverse $\nu$-stable subordinator with respect to $y$, given in [6]. ∎

Lemma 7.2**.**

Let $\sigma_{\nu}(t)$ be a $\nu$-stable subordinator and $T$ be a non-negative random variable independent of $\sigma_{\nu}$. Then

[TABLE]

Proof.

Let us just observe that

[TABLE]

where the second equality is given by $\sigma_{\nu}(t)\overset{d}{=}t^{\frac{1}{\nu}}\sigma_{\nu}(1)$. ∎

By using these two lemmas it is easy to show that

Corollary 7.3**.**

Let $T\sim ML_{\nu}(\lambda)$, $\sigma_{\nu}(t)$ be a $\nu$-stable subordinator and $S\sim Exp(\lambda)$ independent from $\sigma_{\nu}(t)$. Then $T\overset{d}{=}S^{\frac{1}{\nu}}\sigma_{\nu}(1)$.

To simulate the behaviour of the queue process $Q^{\nu}(t)$, we will use a modified version of the Gillespie algorithm (see [13]). This modified version, already introduced in [7], is based on the fact that if $T_{n}$ are the jump times of $Q^{\nu}(t)$, then $Q^{\nu}(T_{n})$ is still a Markov chain with the same transition probabilities as $Q^{1}(T_{n})$. Thus we need first to simulate the jump times $T_{n}$ and then to define what event will happen in such time. In particular we will follow these steps

1. Step 1

Fix $Q^{\nu}(0)$ and $T_{0}=0$; 2. Step 2

Suppose have simulated the queue up to the $n$-th event that happens in time $T_{n}$. If $Q(T_{n})=(0,0)$ simulate a random variable $I\sim ML_{\nu}(\lambda)$, if $Q(T_{n})\in\operatorname{\mathcal{S}}^{*}$ simulate a random variable $I\sim ML_{\nu}(\lambda+k\mu)$; 3. Step 3

Obtain the random variable $T_{n+1}=T_{n}+I$; 4. Step 4

If $Q^{\nu}(T_{n})=(0,0)$, pose $Q^{\nu}(T_{n+1})=(1,k)$. If $Q^{\nu}(T_{n})\in\operatorname{\mathcal{S}}^{*}$, then simulate a uniform variable $U$ in $(0,1)$. If $U<\frac{\lambda}{\lambda+k\mu}$, then pose $N(T_{n+1})=N(T_{n})+1$ and $S(T_{n+1})=S(T_{n})$. If $U\geq\frac{\lambda}{\lambda+k\mu}$, distinguish three cases:

• •

if $S(T_{n})>1$, then pose $N(T_{n+1})=N(T_{n})$ and $S(T_{n+1})=S(T_{n})-1$;

• •

if $S(T_{n})=1$ and $N(T_{n})>1$, then pose $S(T_{n+1})=k$ and $N(T_{n+1})=N(T_{n})-1$;

• •

if $S(T_{n})=1$ and $N(T_{n}))=1$, then pose $S(T_{n+1})=N(T_{n+1})=0$; 5. Step 5

To obtain the whole process $Q^{\nu}(t)$, one has just to pose $Q^{\nu}(t)=Q^{\nu}(T_{n-1})$ for any $t\in[T_{n-1},T_{n})$.

In Figure 1 we provide the simulation of a sample path of the process $N^{\nu}(t)$ for a fractional Erlang queue and we show its embedded Markov chain.

Another way one can use to simulate such process, is to first simulate a classical Erlang queue $Q(t)$, then to simulate an inverse $\nu$-stable subordinator $L_{\nu}(t)$ (some numerical approach to inverse subordinators are described in [35, 36], while some simulation algorithms for the stable subordinator are given in [3]) and finally to combine the two simulated process in $Q^{\nu}(t):=Q(L_{\nu}(t))$ (see, for instance, [29, Example $5.21$]). In this way, the fractional Erlang queue is obtained as a dilatation in time of the classical one, thus we can explicitly see the effect of the inverse $\nu$-stable subordinator, as it is done in Figure 2.

Acknowledgments

N. Leonenko was supported in particular by Australian Research Council’s Discovery Projects funding scheme (project DP160101366), and by project MTM2015-71839-P of MINECO, Spain (co-funded with FEDER funds).

E. Pirozzi was supported by INDAM-GNCS and MANM.