Fractional Erlang Queues
Giacomo Ascione Nikolai Leonenko Enrica Pirozzi

Future frame semantic segmentation of time-lapsed videos with large temporal displacement
     Talha Siddiqui Samarth Bharadwaj

IBM Research Labs, India


    December 30, 2023
=========================================================================================






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.





§ 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 (<cit.>). 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 (<cit.>). In <cit.>, 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 (<cit.>).

The first fractional generalization of the classical M/M/1 queue process was proposed in <cit.>. The transient behaviour of the fractional M/M/1 queues with catastrophes (which in the classical case are studied for instance in <cit.> and further generalized in <cit.>) is investigated in <cit.>.

In this paper, we present a study of the transient behaviour of the fractional Erlang queue M/E_k/1.

In Section <ref>, 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 <ref> 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 <ref> 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 <ref> formulas for transient state probabilities of such process and their Laplace transform are finally obtained.
In Section <ref> some characteristics of the process are studied. In particular in Subsection <ref>, 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 <ref> we provide also a formula for the distribution of the busy period. In Subsection <ref> we focus on the waiting times, determining a formula for a particular conditional waiting time.

Finally, in Section <ref>, we provide an algorithm for sample paths simulation of such process, which is a modified version of the well-known Gillespie algorithm (see <cit.>).


§ 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 λ and an Erlang service system with shape parameter k and rate μ. Such service system can be seen as a system composed of k phases each with exponential service time with parameter kμ.

Let us denote with N(t) the number of the customers in the system at the time t ≥ 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 ≥ 0. If N(t)=0 let us denote also S(t)=0. Finally denote Q(t)=(N(t),S(t)) for t ≥ 0, that is a sort of state-phase process with state space:

    ={(n,s)∈×:  n ≥ 1,  1 ≤ s ≤ k}∪{(0,0)}

and let us denote ^*={(n,s)∈×:  n ≥ 1,  1 ≤ s ≤ k}.
Let us define the transient state-phase probability functions as 

    p_n,s(t)   =(Q(t)=(n,s)|Q(0)=(0,0))     (n,s)∈^*,  t ≥ 0 
    
    p_0(t)   =(Q(t)=(0,0)|Q(0)=(0,0)),  t ≥ 0.

It is well known (see for instance <cit.>) that the function p_n,s(t) and p_0(t) satisfy the following difference-differential equations:

    p_0t(t)=-λ p_0(t)+kμ p_1,1(t) 
    p_1,st(t)=-(λ+kμ)p_1,s(t)+kμ p_1,s+1(t)     1 ≤ s ≤ k-1 
    p_1,kt(t)=-(λ+kμ)p_1,k(t)+kμ p_2,1(t)+λ p_0(t) 
    p_n,st(t)=-(λ+kμ)p_n,s(t)+kμ p_n,s+1(t)+λ p_n-1,s(t)     n≥ 2,  1 ≤ s ≤ k-1 
    p_n,kt(t)=-(λ+kμ)p_n,k(t)+kμ p_n+1,1(t) + λ p_n-1,k(t)     n≥ 2
    
    p_0(0)=1 
    
    p_n,s(0)=0     n ≥ 1,  1 ≤ s ≤ k.

In <cit.> a solution of this system is obtained by means of generalized modified Bessel functions defined in <cit.> as

    I_n^k(t)=(t/2)^n∑_r=0^+∞(t/2)^r(k+1)/r!Γ(n+rk+1),    t ∈,  n=0,1,2,…,  k=1,2,…

and in <cit.> as

    I_n^k,s(t)=(t/2)^n+k-s∑_r=0^+∞(t/2)^r(k+1)/(k(r+1)-s)!Γ(n+r+1)
     t ∈,  n=0,1,2,…,  k=1,2,…,  s=1,2,…,k.

In particular we have:

    p_0(t)=∑_m=1^+∞m(λ/kμ)^-m/k+1I_m^k(2(λ/kμ)^1/k+1kμ t)/kμ te^-(λ+kμ)t

while for n≥ 1 and 1 ≤ s ≤ k-1:

    p_n,s(t)   =(λ/kμ)^k(n-1)+s/k+1I_n^k,s(2(λ(kμ)^k)^1/k+1t)e^-(λ+kμ)t
       +kμ(λ/kμ)^k(n-1)+s/k+1∫_0^tp_0(z)I_n^k,s(2(λ(kμ)^k)^1/k+1(t-z))e^-(λ+kμ)(t-z)dz
       -kμ(λ/kμ)^k(n-1)+s+1/k+1∫_0^tp_0(z)I_n^k,s+1(2(λ(kμ)^k)^1/k+1(t-z))e^-(λ+kμ)(t-z)dz

and finally for n ≥ 1 and s=k:

    p_n,k(t)   =(λ/kμ)^kn/k+1I_n^k,k
    (2(λ(kμ)^k)^1/k+1t)e^-(λ+kμ)t
       +kμ(λ/kμ)^kn/k+1∫_0^tp_0(z)I_n^k,k(2(λ(kμ)^k)^1/k+1(t-z))
    e^-(λ+kμ)(t-z)dz
       -kμ(λ/kμ)^kn+1/k+1∫_0^tp_0(z)
    I_n+1^k,1(2(λ(kμ)^k)^1/k+1(t-z))e^-(λ+kμ)(t-z)dz.

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

    p_0(t)=∑_m=1^+∞∑_r=0^+∞mλ^r(kμ)^m+rk-1/r!Γ(n+rk+1)t^m+r(k+1)-1e^-(λ+kμ)t

while for n ≥ 1 and 1 ≤ s ≤ k-1 we have

    p_n,s(t)   =∑_j=0^+∞λ^n+j(kμ)^k(j+1)-s/(k(j+1)-s)!Γ(n+j+1)t^n+k-s+j(k+1)e^-(λ+kμ)t
       +∑_j=0^+∞λ^n+j(kμ)^k(j+1)-s+1/(k(j+1)-s)!Γ(n+j+1)
       ×∫_0^tp_0(z)(t-z)^n+k-s+j(k+1)e^-(λ+kμ)(t-z)dz
       -∑_j=0^+∞λ^n+j(kμ)^k(j+1)-s/(k(j+1)-s-1)!Γ(n+j+1)
       ×∫_0^tp_0(z)(t-z)^n+k-s-1+j(k+1)e^-(λ+kμ)(t-z)dz

and finally for n≥ 1 and s=k we have

    p_n,k(t)   =∑_j=0^+∞λ^n+j(kμ)^kj/(kj)!Γ(n+j+1)t^n+j(k+1)e^-(λ+kμ)t+
       +∑_j=0^+∞λ^n+j(kμ)^kj+1/(kj)!Γ(n+j+1)×∫_0^tp_0(z)(t-z)^n+j(k+1)e^-(λ+kμ)(t-z)dz
       -∑_j=0^+∞λ^n+j+1(kμ)^k(j+1)/(k(j+1)-1)!
    	Γ(n+j+2)×∫_0^tp_0(z)(t-z)^n+k+j(k+1)e^-(λ+kμ)(t-z)dz.

We can also give an equivalent representation of the queue by using the queue length process, that is to say the process (t) given by:

    (t)= k(N(t)-1)+S(t)     N(t)>0 
    
    0     N(t)=0.

The queue length process gives us the length of the queue by means of phases. In particular, since the map m_k:→_0 given by:

    m_k(n,s)= k(n-1)+s     (n,s)∈^* 
    
    0     (n,s)=(0,0)

is bijective, each state of (t) represents only one state of Q(t) and vice-versa, so that we can equivalently describe the queue with Q(t) or (t): indeed we have m_k(Q(t))=(t). Let us also recall that the inverse map of m_k(n,s) is given by (n_k(m),s_k(m)) where:

    s_k(m)=min{s>0:  s ∈ [m]_k}    m>0 
    
    0     m=0

where [m]_k is the residual class of m modulo k and:

    n_k(m)=m-s_k(m)/k+1     m>0 
    
    0     m=0.

Let us define the state probabilities of (t) as 

    _n(t)=((t)=n|(0)=0).

In <cit.> and <cit.> the forward equations for _n are given as:

    _0t(t)=-λ_0(t)+kμ_1(t) 
    _nt(t)=-(λ+kμ) _n(t)+kμ_n+1(t)+λ∑_m=1^nc_m_n-m(t)     n ≥ 1
    _n(0)=δ_n,0    n ≥ 0

where c_m=δ_m,k for M/E_k/1 queues. However, the bijectivity of m_k(n,s) gives us the equality

    _m(t)=p_n_k(m),s_k(m)(t)

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

In <cit.> the mean queue length, defined as M(t)=[(t)|(0)=0], is show to solve the following Cauchy problem:

    Mt(t)=k(λ-μ)+kμ_0(t)
    
    M(0)=0

which leads, in <cit.>, to the formula:

    M(t)=k(λ-μ)t+kμ∫_0^t_0(y)dy.

Moreover, in <cit.>, the distribution function of the busy period has been obtained as

    B(t)=∑_r=0^+∞kλ^r(kμ)^k(r+1)/r!Γ(rk+k+1)∫_0^t z^k+r(k+1)-1e^-(λ+kμ)zdz.

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 <cit.> as:

    w(ξ;t)=∑_n=0^+∞_n(t)kμ(kμ t)^n-1/(n-1)!e^-kμ t.

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


§ FRACTIONAL ERLANG QUEUES

In <cit.> 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 {σ_ν(t)}_t ≥ 0 a ν-stable subordinator independent from Q(t), such that (see, for instace, <cit.>) 

    [e^-vσ_ν(t)]=e^-tv^ν,  v>0,  ν∈ (0,1).

Define

    L_ν(t)=inf{s ≥ 0:  σ_ν(s)> t}   t≥ 0

the inverse of the ν-stable subordinator and pose N^ν(t):=N(L_ν(t)), S^ν(t):=S(L_ν(t)) and 
    Q^ν(t):=(N^ν(t),S^ν(t))=Q(L_ν(t)).

We will say that the process Q^ν(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^ν(t). Let us define:

    p^ν_n,s(t)   =(Q^ν(t)=(n,s)|Q^ν(0)=(0,0))     (n,s)∈^* 
    
    p^ν_0(t)   =(Q^ν(t)=(0,0)|Q^ν(0)=(0,0)).

To obtain such forward equations, we will need the fractional integral (see <cit.>) of order ν∈ (0,1) defined, for a function x:[0,t_1]⊆→, as

    _t^ν x=1/Γ(ν)∫_0^t(t-τ)^ν-1x(τ)dτ

for any t ∈ [0,t_1] and the Caputo fractional derivative of order ν∈ (0,1)

    D^ν_tx=1/Γ(1-ν)∫_0^txt(τ)dτ/(t-τ)^ν.

The classes of functions for which such operators are well-defined are discussed in <cit.>.

Moreover, we will use a Laplace transform approach, as described in <cit.>. Let us also remark that if x is the Laplace transform of the function x and  is the Laplace transform operator, then:

    [D_t^ν x](z)=z^νx(z)-z^ν-1x(0).

For the state probabilities we can show the following theorem.

The state probabilities (p^ν_n,s)_(n,s)∈ are solution of the following fractional difference-differential Cauchy problem:

    D_t^ν p_0^ν=-λ p_0^ν(t)+kμ p_1,1^ν(t) 
    
    D_t^ν p_1,s^ν=-(λ+kμ)p_1,s^ν(t)+kμ p_1,s+1^ν(t)     1 ≤ s ≤ k-1 
    
    D_t^ν p_1,k^ν=-(λ+kμ)p_1,k^ν(t)+kμ p_2,1^ν(t)+λ p_0^ν(t) 
    
    D_t^ν p_n,s^ν=-(λ+kμ)p_n,s^ν(t)+kμ p_n,s+1^ν(t)+λ p_n-1,s^ν(t)     n ≥ 2,  1 ≤ s ≤ k-1 
    
    D_t^ν p_n,k^ν=-(λ+kμ)p_n,k^ν(t)+kμ p_n+1,1^ν(t)+λ p_n-1,k^ν(t)     n ≥ 2 
    
    p_0^ν(0)=1 
    
    p_n,s^ν(0)=0     n ≥ 1,  1 ≤ s ≤ k.

where D_t^ν is the Caputo fractional derivative.


Let us define the probability generating function

    G^ν(z,t)=∑_n=1^+∞∑_s=1^kz^k(n-1)+sp_n,s^ν(t),  |z|≤ 1.

By multiplying the second equation of (<ref>) by z^s+1, the third equation of (<ref>) by z^k+1, the fourth equation of (<ref>) by z^k(n-1)+s+1, the fifth equation of (<ref>) by z^kn+1 and the summing all these equations with respect to n and s we have that the state probabilities p^ν_n,s solve this Cauchy problem if and only if G^ν(z,t) solves the following Cauchy problem:

    zD_t^ν G^ν(z,·)=(1-z)[G^ν(z,t)(kμ-λ(z+…+z^k))-kμ p_0^ν(t)]
    
    G^ν(z,0)=0.

Denoting with G^ν(z,v) and π_0^ν(v) the Laplace transform of G^ν(z,t) and p_0^ν(t), we know that G^ν(z,t) solves the Cauchy problem (<ref>) if and only if its Laplace transform solves:

    z v^νG^ν(z,v)=(1-z)[G^ν(z,v)(kμ-λ(z+…+z^k))-kμπ_0^ν(v)]

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

    p_n,s^ν(t)=(Q^ν(t)=(n,s)|Q^ν(0)=(0,0))=(Q(L_ν(t))=(n,s)|Q(0)=(0,0))
       =∫_0^+∞(Q(y)=(n,s)|Q(0)==(0,0))(L_ν(t)∈ dy)=∫_0^+∞p_n,s^1(t)(L_ν(t)∈ dy)

and, analogously:

    p_0^ν(t)=∫_0^+∞p_0^1(y)(L_ν(t)∈ dy).

Thus for the probability generating function we know that

    G^ν(z,t)=∫_0^+∞G^1(z,y)(L_ν(t)∈ dy).

Let us recall from <cit.> the following Laplace transform

    [(L_ν(t)∈ dy)](v)=v^ν-1e^-yv^νdy,  v>0.

Thus by using Eq. (<ref>) in Eqs. (<ref>) and (<ref>) we obtain:

    π_0^ν(v)   =∫_0^+∞p_0^1(y)v^ν-1e^-yv^νdy,
    G^ν(z,v)   =∫_0^+∞G^1(z,y)v^ν-1e^-yv^νdy.

Using these two formulas we have that G^ν(z,v) solves Eq. (<ref>) if and only if

    z v^ν   ∫_0^+∞G^1(z,y)e^-yv^νdy=
       =∫_0^+∞(1-z)[G^1(z,y)(kμ-λ(z+…+z^k))-kμ p_0^1(y)]e^-yv^νdy.

From <cit.> we know that G^1(z,t) is solution of the following Cauchy problem:

    zG^1t(z,t)=(1-z)[G^1(z,t)(kμ-λ(z+…+z^k))-kμ p_0^1(t)]
    
    G^1(z,0)=0.

Using such property in (<ref>) we obtain

    zv^ν∫_0^+∞G^1(z,y)e^-yv^νdy=z∫_0^+∞G^1(z,y)te^-yv^νdy

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:

    ^ν(t)= 0     N^ν(t)=0
    
    k(N^ν(t)-1)+S^ν(t)     N^ν(t)≠ 0

and define _m^ν(t):=(^ν(t)=m|^ν(0)=0). By the equality _m^ν(t)=p^ν_n_k(m),s_k(m)(t) and Theorem <ref>, we know that (_m^ν(t))_m ∈ is solution of the fractional Cauchy problem

    D_t^ν_0^ν=-λ_0^ν(t)+kμ_1^ν(t) 
    D_t^ν_n^ν=-(λ+kμ) _n^ν(t)+kμ_n+1^ν(t)+λ∑_m=1^nc_m_n-m^ν(t)     n ≥ 1
    _n^ν(0)=δ_n,0    n ≥ 0

where c_j=δ_j,k, which is equivalent to the fractional Cauchy problem (<ref>). Since the solutions ^ν_m(t) are defined as probabilities, by using Corollary 2 of <cit.> and Schur's test (see <cit.>) one can prove the following:

The fractional Cauchy problems (<ref>) and (<ref>) admit unique global solutions 𝐩^ν(t)=(p_n,s^ν(t))_(n,s) ∈ and 𝐏^ν(t)=(_m^ν(t))_m ∈, where 𝐩^ν,𝐏^ν:[0,+∞) → l^2.




§ INTERPRETATION OF THE FRACTIONAL M/E_K/1 QUEUE


By using the result given in Theorem <ref> 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:

    h^ν_n(t)   :=(N^ν(t)=n|Q^ν(0)=(0,0))  n ≥ 1
    
    q^ν_s(t)   :=(S^ν(t)=s|Q^ν(0)=(0,0))  1 ≤ s ≤ k.

Note that

    h^ν_n(t)=∑_s=1^kp^ν_n,s(t)        q^ν_s(t)=∑_n=1^+∞p^ν_n,s(t).

Let us work with the functions h_n(t). By summing the equations of (<ref>) with respect to s we obtain the following fractional Cauchy problem

    D_t^ν p_0^ν=-λ p_0^ν(t)+kμ p_1,1^ν(t) 
    
    D_t^ν h_1^ν=-λ h_1^ν(t)+kμ(p_2,1^ν(t)-p_1,1^ν(t))+λ p_0^ν(t)
    
    D_t^ν h_n^ν=-λ h_n^ν(t)+kμ(p_n+1,1^ν(t)-p_n,1^ν(t))+λ h_n-1^ν(t)     n ≥ 2
    
    p_0^ν(t)=1 
    
    h_n^ν(t)=0     n ≥ 1

that, for ν=1, gives us

    p_0^1t(t)=-λ p_0^1(t)+kμ p_1,1^1(t) 
    h_1^1t(t)=-λ h_1^1(t)+kμ(p_2,1^1(t)-p_1,1^1(t))+λ p_0^1(t)
    h_n^1t(t)=-λ h_n^1(t)+kμ(p_n+1,1^1(t)-p_n,1^1(t))+λ h_n-1^1(t)     n ≥ 2
    
    p_0^1(t)=1 
    
    h_n^1(t)=0     n ≥ 1.

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

    F_T(t)=1-E_ν(-λ t^ν),  t ≥ 0,  λ >0,

where E_ν is the Mittag-Leffler function (see <cit.>) defined as

    E_ν(z)=∑_k=0^+∞z^k/Γ(ν k+1), ν>0, z ∈;

we call T a Mittag-Leffler random variable with fractional index ν and parameter λ and we denote it by T ∼ ML_ν(λ). In particular let us denote Ψ(t)=E_ν(-λ t^ν), f_T the probability density function of T and 





denoting with f_T its Laplace transform and recalling the Laplace transform of the Mittag-Leffler function that can be obtained from (<ref>) by posing δ=γ=1, we have

    f_T(v)=λ/λ+v^ν,  v >0.

Now we can show the following Theorem.

The interarrival times of a ν-fractional M/E_k/1 queue are i.i.d. Mittag-Leffler random variables with fractional index ν and parameter λ.


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 N^1(t) as a modification of N^1(t) such that the state probabilities 

    h^1_n(t):=(N^1(t)=n)  n ≥ 0

are solution of the following Cauchy problem

    h_m^1t(t)=-λh_m^1(t) 
    h_m+1^1t(t)=λh_m^1
    h_n^1t(t)=0  n≥ 0, n ≠ m,m+1 
    h_n^1(0)=δ_n,m

that is to say that N^1(t) is a birth process (μ=0) that has the same birth rate λ 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:

    h_m^1t(t)=-λh_m^1(t) 
    h_m+1^1t(t)=λh_m^1(t)
    h_m^1(0)=1 
    h_m+1^1(0)=0.

Moreover, since they are probabilities, we have also that h_m^1(t)+h_m+1^1(t)=1 for any t ≥ 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 h_m+1^1(t)=F_T^1(t). 

Now let us define N^ν(t)=N^1(L_ν(t)) and

    h^ν_m(t):=(N^ν(t)=m)       h^ν_m+1(t):=(N^ν(t)=m+1).

By definition, N^ν(t) is a process that starts from m and admits m+1 as absorbent state. Moreover, denoting with T^ν the arrival time of the first new customer in the queue Q^ν starting from N^ν(0)=m, we have that h_m+1^ν(t)=F_T^ν(t). Indeed, since N^1(t) is a Markov process, N^ν(t) is a semi-Markov process and, denoting with T_n the n-th jump time, the discrete time process N^ν(T_n) is a time-homogeneous Markov chain (see, for instance, <cit.>). For this reason, T^ν can be seen as the interarrival time of the queue Q^ν.

Now we have to show that

    D_t^νh_m^ν=-λh_m^ν(t) 
    
    D_t^νh_m+1^ν=λh_m^ν(t)
    h_m^ν(0)=1 
    h_m+1^ν(0)=0.

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 h_m^ν the Laplace transform of h_m^ν:

    v^νh_m^ν(v)-v^ν-1=-λh_m^ν(v),  v>0.

By using the definition of h_m^ν(t) we have that:

    h_m^ν(t)=∫_0^+∞h_m^1(y)(L_ν(t)∈ dy)

and thus, taking the Laplace transform

    h_m^ν(v)=v^ν-1∫_0^+∞h_m^1(y)e^-yv^νdy.

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

    v^ν∫_0^+∞h^1_m(y)e^-yv^νdy-1=∫_0^+∞-λh^1_m(y)e^-yv^νdy.

Recalling that h^1_m(t) is solution of the following Cauchy problem

    h^1_mt(t)=-λh^1_m(t)
    h^1_m(0)=1

we have that Eq. (<ref>) is verified if and only if

    v^ν∫_0^+∞h^1_m(y)e^-yv^νdy-1=∫_0^+∞h^1_m(y)te^-yv^νdy

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

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

    D_t^νh_m^ν=-λh_m^ν(t) 
    h_m^ν(0)=1

whose solution is h_m^ν(t)=E_ν(-λ t^ν) (see <cit.>). Finally:

    F_T^ν(t)=h_m+1^ν(t)=1-h_m^ν(t)=1-E_ν(-λ t^ν),  t ≥ 0.


Now let us work with the functions q_s^ν(t), in order to obtain some information on the service times. Summing the equations in (<ref>) with respect to n we have:

    D_t^ν p_0^ν=-λ p_0^ν(t) + k μ p_1,1^ν(t)
    
    D_t^ν q_s^ν=-kμ q_s^ν(t)+kμ q_s+1^ν(t)     1 ≤ s ≤ k-1 
    
    D_t^ν q_k^ν=-kμ q_k^ν(t)+kμ q_1^ν(t)-kμ p_1,1^ν(t)+λ p_0^ν(t)
    
    p_0^ν(0)=1
    
    q_s^ν(0)=0     1 ≤ s ≤ k

that, for ν=1, gives us

    p_0^1t(t)=-λ p_0^1(t) + k μ p_1,1^1(t)
    q_s^1t(t)=-kμ q_s^1(t)+kμ q_s+1^1(t)     1 ≤ s ≤ k-1 
    q_k^1t(t)=-kμ q_k^1(t)+kμ q_1^1(t)-kμ p_1,1^1(t)+λ p_0^1(t)
    
    p_0^1(0)=1
    
    q_s^1(0)=0     1 ≤ s ≤ k.

We can use the functions q_s^ν exactly in the same way we did with h_n^ν to show the following Theorem.

The inter-phase times of a ν-fractional M/E_k/1 queue are i.i.d. Mittag-Leffler random variables with fractional index ν and parameter kμ.

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

    F_T(t)=1-∑_n=0^k-1(λ t^ν)^n/n!E_ν^(n)(-λ t^ν)=1-∑_n=0^k-1∑_j=1^+∞(-1)^jn+jn(λ t^ν)^j+n/Γ(ν(j+n)+1),  t ≥ 0;

we call T a generalized Erlang random variable (as introduced in <cit.>) with fractional index ν, shape parameter k and rate λ and we denote it by T ∼ GE_ν(k,λ). In particular T∼ GE_ν(k,λ) is sum of k i.i.d. Mittag-Leffler random variable of fractional index ν and parameter λ. For this reason, denoting with f_T the probability density function of T and f_T its Laplace transform, we have, from (<ref>),

    f_T(v)=λ^k/(λ+v^ν)^k,  v>0.

From Theorem <ref> we easily obtain the following Corollary

The service times of a ν-fractional M/E_k/1 queue are i.i.d. generalized Erlang random variables with fractional index ν, shape parameter k and rate kμ.


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 ν and parameter kμ.

From Theorems <ref> and <ref> and Corollary <ref> 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 <cit.>.

We can also obtain the distribution of a general sojourn time of ^ν(t) working with _n^ν(t) as we already did with q_s^ν(t) and h_n^ν(t).

The sojourn time of the queue length process of a ν-fractional M/E_k/1 queue in a state n>0 is a Mittag-Leffler random variable with fractional index ν and parameter λ+kμ.


The sojourn time of ^ν(t) in the state 0 coincides with an interarrival time, so it is a Mittag-Leffler random variable of fractional index ν and parameter λ.

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

Interarrival times and inter-phase times are not independent for ν∈ (0,1).


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

    (min{T,S}≥ t)=(T ≥ t, S ≥ t)=(T ≥ t)(S ≥ t).

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

    (T_S≥ t)=E_ν(-(λ+kμ)t^ν),       (T_P≥ t)=E_ν(-kμ t^ν),       (T_A≥ t)=E_ν(-λ t^ν)

thus we have that (T_S ≥ t)=(T_P ≥ t)(T_A ≥ t) if and only if

    E_ν(-(λ+kμ)t^ν)=E_ν(-kμ t^ν)E_ν(-λ t^ν)

that is true if and only if λ=kμ=0 or ν=1 (see <cit.>).



§ TRANSIENT STATE PROBABILITIES

Our aim is now to obtain the transient state probability functions p_0^ν(t) and p_n,s^ν(t) in a closed form, with the aid of the three-parameter Mittag-Leffler function E_ν,μ^δ (see <cit.>) defined as

    E_ν,μ^δ(z)=∑_k=0^+∞Γ(δ+k)z^k/k!Γ(δ)Γ(ν k+μ),  (γ),(ν),(δ)>0,  z ∈

and the Laplace transform formula (see <cit.>)

    [z^γ-1E_ν,γ^δ(wz^ν)](v)=v^ν,δ-γ/(v^ν-w)^δ, 
    γ,ν,δ,w ∈,  (γ),(ν),(δ)>0,  v ∈, |wv^ν|<1.

To do this, let us first work with p_0^ν(t).

Let p_0^ν(t) be the 0 state probability of a fractional M/E_k/1 queue Q^ν(t) and π_0^ν(v) its Laplace transform. Then we have

    p_0^ν(t)=∑_m=1^+∞∑_r=0^+∞C^0_m,rt^γ^0_m,r-1E_ν,γ^0_m,r^δ^0_m,r(-(λ+kμ)t^ν)

and

    π_0^ν(v)=∑_m=1^+∞∑_r=0^+∞C^0_m,rv^νδ^0_m,r-γ^0_m,r/(λ+kμ+v^ν)^δ^0_m,r

where

    C^0_m,r=m/m+r(k+1)m+r(k+1)m+rkλ^r(kμ)^m+rk-1,
    δ^0_m,r=m+r(k+1),      γ^0_m,r=ν(δ^0_m,r-1)+1.



Let us recall that we already have a closed form for p_0^1(t), given by Eq. (<ref>). Moreover, using this closed form in Eq. (<ref>) we have:

    π_0^ν(v)=∑_m=1^+∞∑_r=0^+∞mλ^r(kμ)^m+rk-1/r!Γ(m+rk+1)v^ν-1∫_0^+∞y^m+r(k+1)-1e^-(λ+kμ+v^ν)ydy

that, integrating, becomes

    π_0^ν(v)=∑_m=1^+∞∑_r=0^+∞mλ^r(kμ)^m+rk-1(m+r(k+1)-1)!/r!Γ(m+rk+1)v^ν-1/(λ+kμ+v^ν)^m+r(k+1).

Now let us pose C^0_m,r, δ^0_m,r and γ^0_m,r as in (<ref>) and observe that

    π_0^ν(v)=∑_m=1^+∞∑_r=0^+∞C^0_m,rv^νδ^0_m,r-γ^0_m,r/(λ+kμ+v^ν)^δ^0_m,r.

Finally, by using Eq. (<ref>) we obtain:

    p_0^ν(t)=∑_r=0^+∞∑_m=1^+∞C^0_m,rt^γ^0_m,r-1E_ν,γ^0_m,r^δ^0_m,r(-(λ+kμ)t^ν).



For k=1 we have that:

    C^0_m,r=m/m+2rm+2rm+rλ^rμ^m+r-1,
    δ^0_m,r=m+2r,     γ^0_m,r=ν(m+2r-1)+1

and then we have

    p_0^ν(t)=∑_r=0^+∞m/m+2rm+2rm+rλ^rμ^m+r-1t^ν(m+2r-1)E^m+2r_ν,ν(m+2r-1)+1(-(λ+μ)t^ν).

Now, let us change variables, posing m=m+r, so we have

    p_0^ν(t)=∑_r=0^+∞∑_m=r+1^+∞m-r/m+rm+rmλ^rμ^m-1t^ν(m+r-1)E^m+r_ν,ν(m+r-1)+1(-(λ+μ)t^ν).

This is the second formula obtained in Remark 4 in <cit.> for the fractional M/M/1 queue. 


For ν=1 we obtain

    γ^0_m,r=δ^0_m,r=m+r(k+1)

and then, with some calculations

    p_0^1(t)   =
    
    ∑_m=1^+∞∑_r=0^+∞C^0_m,r/Γ(δ^0_m,r)t^δ^0_m,r-1e^-(λ+kμ)t.

Let us remark that

    C^0_m,r/Γ(δ^0_m,r)   =mλ^r(kμ)^m+rk-1(m+r(k+1)-1)!/r!Γ(m+rk+1)Γ(m+r(k+1))
    =mλ^r(kμ)^m+rk-1/r!Γ(m+rk+1)

so that

    p_0^1(t)=∑_r=0^+∞mλ^r(kμ)^m+rk-1/r!Γ(m+rk+1)t^m+r(k+1)-1e^-(λ+kμ)t

that is Eq. (<ref>).

To work with p_n,s^ν(t), it could be useful to show the following Lemma.

For any n ∈ we have

    ∫_0^+∞   ∫_0^yp_0^1(z)(y-z)^ne^-(λ+kμ)(y-z)e^-yv^νdzdy
       =∑_m=1^+∞∑_r=0^+∞C^0_m,rn!/(λ+kμ+v^ν)^m+r(k+1)+n+1

where C^0_m,r is defined in (<ref>).


Let us first notice that:

    ∫_0^+∞∫_0^yp_0^1(z)(y-z)^ne^-(λ+kμ)(y-z)e^-yv^νdzdy=
    =∫_0^+∞∫_0^+∞p_0^1(z)(y-z)^ne^-(λ+kμ+v^ν)(y-z)e^-zv^νχ_[0,y](z)dzdy

and then we can use Fubini's theorem to obtain

    ∫_0^+∞∫_0^+∞p_0^1(z)(y-z)^ne^-(λ+kμ+v^ν)(y-z)e^-zv^νχ_[0,y](z)dzdy=
    =∫_0^+∞∫_0^+∞p_0^1(z)(y-z)^ne^-(λ+kμ+v^ν)(y-z)e^-zv^νχ_[0,y](z)dydz.

Now let us pose w=y-z to obtain

    ∫_0^+∞∫_0^+∞p_0^1(z)(y-z)^ne^-(λ+kμ+v^ν)(y-z)e^-zv^νχ_[0,y](z)dydz
    =∫_0^+∞∫_-z^+∞p_0^1(z)w^ne^-(λ+kμ+v^ν)we^-zv^νχ_[0,z+w](z)dwdz

that, since χ_[0,z+w](z)=χ_[0,+∞)(w), gives us

    ∫_0^+∞∫_-z^+∞p_0^1(z)w^ne^-(λ+kμ+v^ν)we^-zv^νχ_[0,z+w](z)dwdz=
    =(∫_0^+∞p_0^1(z)e^-zv^νdz)(∫_0^+∞w^ne^-(λ+kμ+v^ν)wdw).

By integrating, we have

    ∫_0^+∞w^ne^-(λ+kμ+v^ν)wdw=n!/(λ+kμ+v^ν)^n+1.

Moreover, by (<ref>) and (<ref>) we have

    ∫_0^+∞   p_0^1(z)e^-zv^νdz=π_0^ν(v)/v^ν-1
    =∑_m=1^+∞∑_r=0^+∞C^0_m,r/(λ+kμ+v^ν)^m+r(k+1)

so finally we obtain

    ∫_0^+∞   ∫_0^yp_0^1(z)(y-z)^ne^-(λ+kμ)(y-z)e^-yv^νdzdy
       =∑_m=1^+∞∑_r=0^+∞C^0_m,rn!/(λ+kμ+v^ν)^m+r(k+1)+n+1.


Now we are ready to show the following Theorem.

For any (n,s)∈^* let p_n,s^ν(t)=(Q^ν(t)=(n,s)|Q^ν(0)=(0,0)) and π_n,s^ν(v) be its Laplace transform. Then

    p_n,s^ν(t)   =∑_j=0^+∞A^n,s_jt^α_j^n,s-1E^a^n,s_j_ν,α_j^n,s(-(λ+kμ)t^ν)
       +∑_j=0^+∞∑_m=1^+∞∑_r=0^+∞B^n,s_j,m,rt^β_j,m,r^n,s-1E^b^n,s_j,m,r_ν,β_j,m,r^n,s(-(λ+kμ)t^ν)
       -∑_j=0^+∞∑_m=1^+∞∑_r=0^+∞C^n,s_j,m,rt^γ_j,m,r^n,s-1E^c^n,s_j,m,r_ν,γ_j,m,r^n,s(-(λ+kμ)t^ν)

and

    π_n,s^ν(v)   =∑_j=0^+∞A^n,s_jv^ν a^n,s_j-α^n,s_j/(λ+kμ+v^ν)^a^n,s_j +∑_j=0^+∞∑_m=1^+∞∑_r=0^+∞B^n,s_j,m,rv^ν b^n,s_j,m,r-β^n,s_j,m,r/(λ+kμ+v^ν)^b^n,s_j,m,r
       -∑_j=0^+∞∑_m=1^+∞∑_r=0^+∞C^n,s_j,m,rv^ν c^n,s_j,m,r-γ^n,s_j,m,r/(λ+kμ+v^ν)^c^n,s_j,m,r

where

    A^n,s_j=n+k+kj+j-sn+jλ^n+j(kμ)^k(j+1)-s
    
    a^n,s_j=n-s+(j+1)(k+1)     α^n,s_j=ν(a^n,s_j-1)+1
    
    B^n,s_j,m,r=kμ C^0_m,rA^n,s_j      b^n,s_j,m,r=δ^0_m,r+a^n,s_j     β^n,s_j,m,r=ν(b^n,s_j,m,r-1)+1
    
    C^n,s_j,m,r=
    kμ C^0_m,rA^n,s+1_j     s≠ k
    
    kμ C^0_m,rA^n+1,1_j     s= k
         c^n,s_j,m,r=δ^0_m,r+a^n,s+1_j     s ≠ k
    δ^0_m,r+a^n+1,1_j     s=k
    
    γ^n,s_j,m,r=ν(c^n,s_j,m,r-1)+1.



From Eq. (<ref>) we easily obtain

    π_n,s^ν(v)=v^ν-1∫_0^+∞p^1_n,s(y)e^-yv^νdy,  v>0

Let us first consider 1 ≤ s ≤ k-1; by using Eq. (<ref>), we have

    π_n,s^ν(v)   =∑_j=0^+∞λ^n+j(kμ)^k(j+1)-s/(k(j+1)-s)!Γ(n+j+1)v^ν-1∫_0^+∞y^n+k-s+j(k+1)e^-(λ+kμ+v^ν)ydy
       +∑_j=0^+∞λ^n+j(kμ)^k(j+1)-s+1/(k(j+1)-s)!Γ(n+j+1)v^ν-1
       ×∫_0^+∞∫_0^yp_0^1(z)(y-z)^n+k-s+j(k+1)e^-(λ+kμ)(y-z)e^-yv^νdzdy
       -∑_j=0^+∞λ^n+j(kμ)^k(j+1)-s/(k(j+1)-s-1)!Γ(n+j+1)v^ν-1
       ×∫_0^+∞∫_0^yp_0^1(z)(y-z)^n+k-s-1+j(k+1)e^-(λ+kμ)(y-z)e^-yv^νdzdy.

Integrating, using Eq. (<ref>) and posing A^n,s_j, a^n,s_j, α^n,s_j, B^n,s_j,m,r, b^n,s_j,m,r, β^n,s_j,m,r, C^n,s_j,m,r, c^n,s_j,m,r and γ^n,s_j,m,r as in equation (<ref>) for s ≠ k we obtain

    π_n,s^ν(v)   =∑_j=0^+∞A^n,s_jv^ν a^n,s_j-α^n,s_j/(λ+kμ+v^ν)^a^n,s_j +∑_j=0^+∞∑_m=1^+∞∑_r=0^+∞B^n,s_j,m,rv^ν b^n,s_j,m,r-β^n,s_j,m,r/(λ+kμ+v^ν)^b^n,s_j,m,r
       -∑_j=0^+∞∑_m=1^+∞∑_r=0^+∞C^n,s_j,m,rv^ν c^n,s_j,m,r-γ^n,s_j,m,r/(λ+kμ+v^ν)^c^n,s_j,m,r.

Finally, by using Eq. (<ref>) we obtain

    p_n,s^ν(t)   =∑_j=0^+∞A^n,s_jt^α_j^n,s-1E^a^n,s_j_ν,α_j^n,s(-(λ+kμ)t^ν)
       +∑_j=0^+∞∑_m=1^+∞∑_r=0^+∞B^n,s_j,m,rt^β_j,m,r^n,s-1E^b^n,s_j,m,r_ν,β_j,m,r^n,s(-(λ+kμ)t^ν)
       -∑_j=0^+∞∑_m=1^+∞∑_r=0^+∞C^n,s_j,m,rt^γ_j,m,r^n,s-1E^c^n,s_j,m,r_ν,γ_j,m,r^n,s(-(λ+kμ)t^ν).

Now, for s=k, from Eq. (<ref>) and (<ref>) we obtain

    π_n,k^ν(v)   =∑_j=0^+∞λ^n+j(kμ)^kj/(kj)!Γ(n+j+1)v^ν-1∫_0^+∞y^n+j(k+1)e^-(λ+kμ+v^ν)ydy
       +∑_j=0^+∞λ^n+j(kμ)^kj+1/(kj)!Γ(n+j+1)v^ν-1
       ×∫_0^+∞∫_0^yp_0^1(z)(y-z)^n+j(k+1)e^-(λ+kμ)(y-z)e^-yv^νdzdy
       -∑_j=0^+∞λ^n+j+1(kμ)^k(j+1)/(k(j+1)-1)!Γ(n+j+2)v^ν-1
       ×∫_0^+∞∫_0^yp_0^1(z)(y-z)^n+k+j(k+1)e^-(λ+kμ)(y-z)e^-yv^νdzdy.

Integrating, using Eq. (<ref>) and posing A^n,s_j, a^n,s_j, α^n,s_j, B^n,s_j,m,r, b^n,s_j,m,r, β^n,s_j,m,r, C^n,s_j,m,r, c^n,s_j,m,r and γ^n,s_j,m,r as in equation (<ref>) for s = k we obtain

    π_n,k^ν(v)   =∑_j=0^+∞A^n,k_jv^ν a^n,k_j-α^n,k_j/(λ+kμ+v^ν)^a^n,k_j
       +∑_j=0^+∞∑_m=1^+∞∑_r=0^+∞B^n,k_j,m,rv^ν b^n,k_j,m,r-β^n,k_j,m,r/(λ+kμ+v^ν)^b^n,k_j,m,r-∑_j=0^+∞∑_m=1^+∞∑_r=0^+∞C^n,k_j,m,rv^ν c^n,k_j,m,r-γ^n,k_j,m,r/(λ+kμ+v^ν)^c^n,k_j,m,r.

Finally, by using Eq. (<ref>) we obtain

    p_n,k^ν(t)   =∑_j=0^+∞A^n,k_jt^α_j^n,k-1E^a^n,k_j_ν,α_j^n,k(-(λ+kμ)t^ν)
       +∑_j=0^+∞∑_m=1^+∞∑_r=0^+∞B^n,k_j,m,rt^β_j,m,r^n,k-1E^b^n,k_j,m,r_ν,β_j,m,r^n,k(-(λ+kμ)t^ν)
       -∑_j=0^+∞∑_m=1^+∞∑_r=0^+∞C^n,k_j,m,rt^γ_j,m,r^n,k-1E^c^n,k_j,m,r_ν,γ_j,m,r^n,k(-(λ+kμ)t^ν).
























§ 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.


 §.§ Mean queue length

To determine the mean queue length, let us consider the queue length process ^ν(t) and define:

    M^ν(t)=[^ν(t)|^ν(0)=0].

Let us first show the following Theorem

The function M^ν(t) is solution of the following fractional Cauchy problem:

    D_t^ν M^ν=k(λ-μ)+kμ^ν_0(t) 
    
    M^ν(0)=0.



The initial condition follows obviously from the definition of M^ν(t).

To obtain the equation, recall that, by definition:

    M^ν(t)=∑_n=1^+∞n_n^ν(t).

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

    D_t^ν M^ν(t)=-(λ+kμ)M^ν(t)+kμ∑_n=1^+∞n_n+1^ν(t)+λ∑_n=1^+∞∑_m=1^nc_m^ν_n-m(t).

Let us observe that

    ∑_n=1^+∞n_n+1^ν(t)=M^ν(t)-∑_n=0^+∞_n+1^ν(t)=M^ν(t)+_0^ν(t)-1

and

    ∑_n=1^+∞∑_m=1^nnc_m_n-m^ν(t)=∑_m=1^+∞∑_n=m^+∞nc_m_n-m^ν(t).

Now let us pose r=n-m to obtain

    ∑_m=1^+∞∑_n=m^+∞nc_m_m-n^ν(t)   =∑_m=1^+∞∑_r=0^+∞(r+m)c_m_r^ν(t)
       =∑_m=1^+∞c_m∑_r=0^+∞r_r^ν(t)+∑_m=1^+∞mc_m∑_r=0^+∞_r^ν(t).

Recalling that c_m=δ_m,k, we have

    ∑_m=1^+∞c_m∑_r=0^+∞r_r^ν(t)+∑_m=1^+∞mc_m∑_r=0^+∞_r^ν(t)=M^ν(t)+k.

Using Eqs. (<ref>) and (<ref>) in (<ref>) we obtain the first equation of (<ref>).

By integrating (<ref>), we obtain the following representation of the mean queue length:

    M^ν(t)=k(λ-μ)t^ν/Γ(ν+1)+kμ^ν_t _0^ν

For k=1 we obtain the mean queue length of a fractional M/M/1 queue, as stated in Theorem 2.3 of <cit.>. For ν=1 we obtain Eq. (<ref>).

Starting from Eq. (<ref>) we can also obtain a closed form for M^ν(t).

Let M^ν(t)=[^ν(t)|^ν(0)=0] and denote with M^ν(v) its Laplace transform. Then

    M^ν(t)=k(λ-μ)t^ν/Γ(ν+1)+kμ∑_m=1^+∞∑_r=0^+∞C^0_m,rt^γ^M_m,r-1E^δ^0_m,r_ν,γ^M_m,r(-(λ+kμ)t^ν)

and

    M^ν(v)=k(λ-μ)1/v^ν+1+kμ∑_m=1^+∞∑_r=0^+∞C^0_m,rv^νδ^0_m,r-γ^M_m,r/(λ+kμ+v^ν)^δ^0_m,r

where

    γ^M_m,r=νδ^0_m,r+1

and C^0_m,r and δ^0_m,r are defined in (<ref>).


From Eq. (<ref>) and the equality _m^ν=p_n_k(m),s_k(m)^ν we easily obtain that

    M^ν(t)=∫_0^+∞M^1(y)(L_ν(t)∈ dy).

Using (<ref>) we obtain

    M^ν(t)=k(λ-μ)∫_0^+∞y(L_ν(t)∈ dy)+kμ∫_0^+∞∫_0^y_0^1(z)dz(L_ν(t)∈ dy).

Taking the Laplace transform we obtain

    M^ν(v)=k(λ-μ)v^ν-1∫_0^+∞ye^-yv^νdy+kμ v^ν-1∫_0^+∞∫_0^y^1_0(z)e^-yv^νdzdy.

For the first integral we easily obtain

    ∫_0^+∞ye^-yv^νdy=1/v^2ν.

For the second one we have

    ∫_0^+∞∫_0^y^1_0(z)e^-yv^νdzdy   =∫_0^+∞∫_z^+∞^1_0(z)e^-yv^νdydz
       =1/v^ν∫_0^+∞_0^1(z)e^-zv^νdz.

By using Eq. (<ref>) we obtain

    ∫_0^+∞_0^1(z)e^-zv^νdz   =∑_m=1^+∞∑_r=0^+∞mλ^r(kμ)^m+rk-1/r!Γ(n+rk+1)∫_0^+∞z^m+r(k+1)-1e^-(λ+kμ+v^ν)zdz
       =∑_m=1^+∞∑_r=0^+∞C^0_m,r1/(λ+kμ+v^ν)^m+r(k+1).

Using Eqs. (<ref>), (<ref>) and (<ref>) in (<ref>) we obtain

    M^ν(v)=k(λ-μ)1/v^ν+1+kμ∑_m=1^+∞∑_r=0^+∞C^0_m,rv^-1/(λ+kμ+v^ν)^m+r(k+1).

Posing then γ^M_m,r as in (<ref>) we obtain

    M^ν(v)=k(λ-μ)1/v^ν+1+kμ∑_m=1^+∞∑_r=0^+∞C^0_m,rv^νδ^0_m,r-γ^M_m,r/(λ+kμ+v^ν)^δ^0_m,r.

Finally, taking the inverse Laplace transform and recalling eq. (<ref>) we obtain

    M^ν(t)=k(λ-μ)t^ν/Γ(ν+1)+kμ∑_m=1^+∞∑_r=0^+∞C^0_m,rt^γ^M_m,r-1E^δ^0_m,r_ν,γ^M_m,r(-(λ+kμ)t^ν).


Comparing Eq. (<ref>) with (<ref>) we also obtain some information on the fractional integral of _0^ν(t). In particular

    ^ν_t _0^ν=∑_m=1^+∞∑_r=0^+∞C^0_m,rt^γ^M_m,r-1E^δ^0_m,r_ν,γ^M_m,r(-(λ+kμ)t^ν).



 §.§ 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 <cit.>.

Let B^ν be the duration of the busy period of a fractional M/E_k/1 queue, B^ν(t) be its distribution function and B^ν(v) be its Laplace transform. Then we have

    B^ν(t)=∑_r=0^+∞C^B_rt^γ^B_r-1E^δ^B_r_ν,γ^B_r(-(λ+kμ)t^ν)

and

    B^ν(v)=∑_r=0^+∞C^B_rv^νδ^B_r-γ^B_r/(λ+kμ+v^ν)^δ^B_r

where

    C^B_r=kμ C^0_k,r      γ^B_r=νδ^0_k,r+1

with C^0_k,r and δ^0_k,r defined in Eqs. (<ref>). 


Let us consider the process ^1(t) whose state probabilities

    _n^1(t)=(^1(t)=n|^1(0)=1)

satisfy the following Cauchy problem:

    _0t(t)=kμ_1(t) 
    _1t(t)=-(λ+kμ) _1(t)+kμ_2(t) 
    _nt(t)=-(λ+kμ) _n(t)+kμ_n+1(t)+λ∑_m=1^nc_m_n-m(t)     n ≥ 2
    _n(0)=δ_n,k    n ≥ 0.

The process ^1(t)behaves as ^1(t) except that it starts from k (that is to say when the first customer entered the queue) instead that from 0 and admits 0 as absorbent state. In particular, by construction, _0^1(t)=B^1(t).

Now let us consider ^ν(t):=^1(L_ν(t)). As done with the interarrival and the inter-phase times, if T_n is the n-th jump time of the process ^ν(t), then ^ν(T_n) is a time-homogeneous Markov chain, then, by construction, _0^ν(t)=B^ν(t).

From the relation between ^1 and ^ν, we have that

    B^ν(t)=∫_0^+∞B^1(y)(L_ν(t)∈ dy).

Thus, by using Eq. (<ref>) we have

    B^ν(t)=∑_r=0^+∞kλ^r(kμ)^k(r+1)/r!Γ(rk+k+1)∫_0^+∞∫_0^yz^k+r(k+1)-1e^-(λ+kμ)zdz(L_ν(t)∈ dy).

Taking the Laplace transform we have

    B^ν(v)=∑_r=0^+∞kλ^r(kμ)^k(r+1)/r!Γ(rk+k+1)v^ν-1∫_0^+∞∫_0^yz^k+r(k+1)-1e^-(λ+kμ)ze^-yv^νdzdy.

Now let us consider the integral and observe that:

    ∫_0^+∞∫_0^yz^k+r(k+1)-1   e^-(λ+kμ)ze^-yv^νdzdy
       =∫_0^+∞∫_z^+∞z^k+r(k+1)-1e^-(λ+kμ)ze^-yv^νdydz
       =1/v^ν∫_0^+∞z^k+r(k+1)-1e^-(λ+kμ+v^ν)zdz
    =1/v^ν(k+r(k+1)-1)!/(λ+kμ+v^ν)^k+r(k+1).

Using Eq. (<ref>) in Eq. (<ref>) we obtain

    B^ν(v)=∑_r=0^+∞kλ^r(kμ)^k(r+1)(k+r(k+1)-1)!/r!Γ(rk+k+1)v^-1/(λ+kμ+v^ν)^k+r(k+1).

Thus, posing C^B_r and γ^B_r as in (<ref>) we obtain

    B^ν(v)=∑_r=0^+∞C^B_rv^νδ^B_r-γ^B_r/(λ+kμ+v^ν)^δ^B_r.

Finally, taking the inverse Laplace transform and using Eq. (<ref>), we obtain

    B^ν(t)=∑_r=0^+∞C^B_rt^γ^B_r-1E^δ^B_r_ν,γ^B_r(-(λ+kμ)t^ν).




 §.§ Distribution of some conditioned waiting times

For classical single-channel queues, a method to obtain the waiting time distribution has been presented in <cit.>. 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 ∼ ML_ν(λ) and fix t_0>0. Observe that

    (T≤ t_0+t|T ≥ t_0)=1-E_ν(-λ(t_0+t)^ν)/E_ν(-λ t_0^ν).

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, <cit.>). This observation leads us to the definition of a new distribution.

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

    F_T(t)=1-E_ν(-λ(t_0+t)^ν)/E_ν(-λ t_0^ν);

we call such random variable residual Mittag-Leffler random variable starting from t_0 with fractional index ν and parameter λ and we denote it by T ∼ RML_ν(t_0,λ). Denote

    Φ(t):=E_ν(-λ(t_0+t)^ν)

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

    Φ(t)=∑_k=0^+∞(-λ)^k/Γ(ν k+1)(t_0+t)^ν k.

Denote with Φ(v) its Laplace transform and observe that

    Φ(v)=∑_k=0^+∞(-λ)^k/Γ(ν k+1)e^vt_0v^-1-ν kΓ(ν k+1,vt_0),  v>0,

where Γ(z,s) is the analytic extension of the upper incomplete gamma function on ^2 (see, for instance, <cit.>). Denote with f_T(t) the probability density function of T and observe that

    f_T(t)=-1/E_ν(-λ t_0^ν)Φ(t)t.

Denoting with f_T(v) its Laplace transform, we obtain

    f_T(v)=1-e^vt_0/E_ν(-λ t_0^ν)∑_k=0^+∞(-λ)^k/Γ(ν k+1)v^-ν kΓ(ν k+1,vt_0),  v>0.

Finally let us also observe that if T ∼ RML_ν(t_0,λ) for ν=1, then T ∼ Exp(λ). Now we are ready to determine such waiting times.

Let W^ν_t be the waiting time of a customer that enters the queue at time t>0. Let us define the random set

    _t^ν(ω)={0 ≤ s ≤ t:  ^ν(s-)(ω)=^ν(s)(ω)+1}

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

    _t^ν=|^ν_t|

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

    Θ_t^ν=sup^ν_t    _t^ν>0
    
    t    _t^ν=0

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

    w^ν(ξ;t,t_0,n)dξ=(W^ν_t ∈ dξ| Θ_t^ν=t_0, ^ν(t)=n).

Since we have conditioned this probability with ^ν(t)=n, we know that a customer has to wait n independent phases to be terminated. Thus we have

    [W^ν_t|Θ_t^ν=t_0, ^ν(t)=n]=∑_j=1^nT_j

where T_j are n independent random variables. For 1 ≤ j ≤ n-1, these are just Mittag-Leffler random variables with T_j ∼ ML_ν(kμ), thus, denoting W=∑_j=1^n-1T_j, we know that W∼ GE_ν(n-1,kμ). For T_n, we know that the last time a phase has been completed is t_0, since we have conditioned the probability with Θ_t^ν=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 ∼ RML_ν(t-t_0,kμ). Thus we have that

    w^ν(ξ;t,t_0,n)=f_W∗ f_T_n(ξ).

In particular we can obtain the Laplace transform w^ν(v;t,t_0,n) of w^ν(ξ;t,t_0,n) as

    w^ν(v;t,t_0,n)=(kμ)^n-1/(kμ+v^ν)^n-1×
    ×[1-e^v(t-t_0)/E_ν(-λ (t-t_0)^ν)∑_k=0^+∞(-kμ)^k/Γ(ν k+1)v^-ν kΓ(ν k+1,v(t-t_0))].

For ν=1, we know that T_n ∼ Exp(kμ) and W∼ E_n-1(kμ), thus T_n+W∼ E_n(kμ) and we know that

    w^1(ξ;t,t_0,n)=(kμ)^nξ^n-1e^-kμξ/(n-1)!=w^1(ξ;t,n)

since the density is independent from t_0. Thus we obtain

    w^1(ξ;t)=∑_n=0^+∞_n^1(t)w^1(ξ;t,n)

that is Eq. (<ref>), that is the one given in <cit.>.


§ 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.

Let T ∼ Exp(λ) and σ_ν(t) be a ν-stable subordinator independent of T. Then σ_ν(T)∼ ML_ν(λ)


Let us just observe that

    (σ_ν(T)>t)   =(T>L_ν(t))
       =∫_0^+∞(T>y)(L_ν(t)∈ dy)
       =∫_0^+∞e^-λ y(L_ν(t)∈ dy)
       =E_ν(-λ t^ν)

where the last equality is given by the Laplace transform of the density (L_ν(t)∈ dy) of the inverse ν-stable subordinator with respect to y, given in <cit.>.


Let σ_ν(t) be a ν-stable subordinator and T be a non-negative random variable independent of σ_ν. Then

    σ_ν(T)d=T^1/νσ_ν(1)



Let us just observe that

    (σ_ν(T)>t)   =∫_0^+∞(σ_ν(y)>t)(T ∈ dy)
       =∫_0^+∞(y^1/νσ_ν(1)>t)(T ∈ dy)
       =(T^1/νσ_ν(1))

where the second equality is given by σ_ν(t)d=t^1/νσ_ν(1).

By using these two lemmas it is easy to show that

Let T ∼ ML_ν(λ), σ_ν(t) be a ν-stable subordinator and S ∼ Exp(λ) independent from σ_ν(t). Then T d=S^1/νσ_ν(1).

To simulate the behaviour of the queue process Q^ν(t), we will use a modified version of the Gillespie algorithm (see <cit.>). This modified version, already introduced in <cit.>, is based on the fact that if T_n are the jump times of Q^ν(t), then Q^ν(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


  Step 1 Fix Q^ν(0) and T_0=0;

  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 ∼ ML_ν(λ), if Q(T_n) ∈^* simulate a random variable I ∼ ML_ν(λ+kμ);

  Step 3 Obtain the random variable T_n+1=T_n+I;

  Step 4 If Q^ν(T_n)=(0,0), pose Q^ν(T_n+1)=(1,k). If Q^ν(T_n) ∈^*, then simulate a uniform variable U in (0,1). If U<λ/λ+kμ, then pose N(T_n+1)=N(T_n)+1 and S(T_n+1)=S(T_n). If U≥λ/λ+kμ, 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;


  Step 5 To obtain the whole process Q^ν(t), one has just to pose Q^ν(t)=Q^ν(T_n-1) for any t ∈ [T_n-1,T_n).

In Figure <ref> we provide the simulation of a sample path of the process N^ν(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 ν-stable subordinator L_ν(t) (some numerical approach to inverse subordinators are described in <cit.>, while some simulation algorithms for the stable subordinator are given in <cit.>) and finally to combine the two simulated process in Q^ν(t):=Q(L_ν(t)) (see, for instance, <cit.>). 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 ν-stable subordinator, as it is done in Figure <ref>. 



§ 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.
plain


