Fractional Queues with Catastrophes and their Transient Behaviour

Fractional Queues with Catastrophes and their Transient Behaviour Giacomo Ascione 1, Nikolai Leonenko 2 and Enrica Pirozzi 1,∗ 1 Università degli Studi di Napoli Federico II, Dipartimento di Matematica e Applicazioni "Renato Caccioppoli", Italy; giacomo.ascione@unina.it, enrica.pirozzi@unina.it 2 Cardiff University, United Kingdom; LeonenkoN@cardiff.ac.uk * Correspondence: enrica.pirozzi@unina.it; Academic Editor: name Version August 30, 2018 submitted to Mathematics Abstract: Starting from the definition of fractional M/M/1 queue given in [3] and M/M/1 queue 1 with catastrophes given in [6], we define and study a fractional M/M/1 queue with catastrophes. 2 In particular we focus our attention on the transient behaviour, in which the time-change plays a 3 key role. We first specify the conditions for the global uniqueness of solutions of the corresponding 4 linear fractional differential problem. Then we provide an alternative expression for the transient 5 distribution of the fractional M/M/1 model, the state probabilities for the fractional queue with 6 catastrophes, the distributions of the busy period for fractional queues without and with catastrophes 7 and, finally, the distribution of the time of the first occurrence of a catastrophe. 8


Introduction
Stochastic models for queueing systems have a wide range of applications in computer systems, sales points, telephone or telematic systems and also in several areas of science including biology, medicine and many others.The well known M/M/1 queueing model [1][2][3][4][5] constitutes the theoretical basis for building many other refined models for service systems.
Due to the Markov nature of its inter-arrival times of the customers and of its service times, the model can be mathematically treated in a simple manner, and, for this reason, it is widely used in many modeling contexts.Nevertheless, in the past few decades, the advent of fractional operators, such as fractional derivatives and integrals (see, for instance, [6] and [7] and references therein), has made it clear that different time scales, themselves random, that preserve memory (therefore not Markovian), allow the construction of more realistic stochastic models.
The introduction of the fractional Caputo derivative into the system of differential-difference equations for an M/M/1-type queue was done in [8], where, for a fractional M/M/1 queue, the state probabilities were determined.In this kind of queue model, the inter-arrival times and service times are characterized by Mittag-Leffler distributions [9]; in this case, the model does not have the property of memory loss that is typical of the exponential distributed times of the classical M/M/1 model.Indeed, a time-changed birth-death process [10,11], by means of an inverse stable subordinator [12], solves the corresponding fractional system of differential-difference equations and fractional Poisson processes [13] characterize the inter-arrival and service times.
The fractional M/M/1 model in [8] is an interesting and powerful model, not only because it is a generalization of the classical one, where the fractional order is set to 1, but also because its range of applications is extremely wide.Its importance can be further augmented by including in the model the occurrence of catastrophes, as it was considered in [14] for the classical M/M/1.
The catastrophe is a particular event that occurs in a random time leading to the instantaneous emptying of the system, or to a momentary inactivation of the system, as, for example, the action of a virus program that can make a computer system inactive [15]; other applications of models with catastrophes can be found in population dynamics and reliability contexts (see [16] and references therein).
Motivated by the mathematical need to enrich the fractional M/M/1 model of [8] with the inclusion of catastrophes, we study in this paper the above model; specifically, we determine the transient distribution, the distribution of the busy period (including that of the fractional M/M/1 queue of [8]) and the probability distribution of the time of the first occurrence of the catastrophe.
For these purposes, we need to guarantee the global uniqueness of the solution of the considered linear fractional Cauchy problem on Banach spaces.After recalling the definitions and known results in Section 2, we address the problem of uniqueness in Section 3. In Section 2, we also provide the transient distribution of the fractional M/M/1 model in an alternative form to that given in [8].In Section 4, the distribution of the busy period for the fractional M/M/1 queue (without catastrophes) is obtained.Here, the time-changed birth-death process plays a key role to derive the results.In Section 5, we define the fractional queue with catastrophes; we are able to obtain the distribution of the transient state probabilities by following a strategy similar to that in [14].We also found the distribution of the busy period and of the time of the first occurrence of the catastrophe starting from the empty system.Some special operators and functions used in this paper are specified in the Appendices A and B.

Definition of a Fractional Process Related to M/M/1 Queues
The classical M/M/1 queue process N(t), t ≥ 0 can be described as continuous time Markov chain whose state space is the set {0, 1, 2, . . .} and the state probabilities satisfy the following differential-difference equations: where δ n,0 is the Kroeneker delta symbol, D t = d dt and α, β > 0 are the entrance and service rates, respectively.
Consider the inverse ν-stable subordinator For 0 < ν < 1, the fractional M/M/1 queue process N ν (t), t ≥ 0 is defined by a non-Markovian time change L ν (t) independent of N(t), t ≥ 0, i.e., This process was defined in [8] and it is non-Markovian with non-stationary and non-independent increments.For ν = 1, by definition, N 1 (t) = N(t), t ≥ 0.Then, for a fixed ν ∈ (0, 1], the state probabilities of the number of customers in the system at time t in the fractional M/M/1 queue are characterized by arrivals and services determined by fractional Poisson processes of order ν ∈ (0, 1] [13] with parameters α and β.They are solutions of the following system of differential-difference equations where C 0 D ν t is the Caputo fractional derivative (see Appendix A).
Using Equation ( 5) and representation (3), the state probabilities are obtained in [8]: as well as its Laplace transform In Equation ( 6), the functions E ρ ν,µ are generalized Mittag-Leffler functions (see Appendix B).Note that p ν n (t) ≥ 0 ∀n ≥ 0 and x ≥ 0, be the density of L ν (t); then it is known (see, i.e., [17]) that and (see, i.e., [18], Proposition 4.1) Using (7) and an analytical expression for p ν n (t) given in [19], we can write down an alternative expression for (6) as where , and Actually, it is easy to see from (7) that thus, using [19] and (3), we have and formula (9) follows.On the other hand, using (8), we have where

Linear Fractional Cauchy Problems on Banach Spaces
In order to describe the transient probabilities for our queues, we will need some uniqueness results for solutions of linear fractional Cauchy problems defined on Banach spaces.To do that, let us recall the following Theorem (Theorem 3.19 from [20]): Theorem 1.Let (X, | • |) be a Banach space and J = [0, T] for some T > 0. Consider the ball B R = {x ∈ X : |x| ≤ R}.Let ν ∈ (0, 1) and f : J × B R → X and consider the following Cauchy problem: where C 0 D ν t is the Caputo derivative operator (see Appendix A).Suppose that: for all x ∈ B R and t ∈ J and such that
The previous theorem can be easily adapted to the case in which J = [t 0 , T + t 0 ] and the starting point of the derivative is t 0 .Since we are interested in linear (eventually non-homogeneous) equations, let us show how the previous theorem can be adapted in such a case.Corollary 1.Consider the system (10) and suppose f (t, x) = Ax + ξ where A : X → X is a linear and continuous operator and ξ ∈ X.Then, there exists a R > |x 0 | and T > 0 such that the system admits a unique solution x * ∈ C ν (J, B R ).
Let us choose T such that the conditions of Theorem 1 are verified.To do that, consider and observe that Thus, one can fix L = 2M( R) Γ(ν+1) and L 0 = A .Moreover, since for fixed ν 1 ∈ (0, ν) the function τ → L A (τ) is decreasing and lim τ→0 L A (τ) = 0, then one can easily find a τ > 0 such that L A (τ) < 1.Since we are under the hypotheses of Theorem 1, then we have shown the local existence and uniqueness of a solution x * ∈ C ν (J, B R).
However, using such corollary, we can only afford local uniqueness.Global uniqueness of the solution of the Cauchy problem (10) can be obtained with the additional hypothesis that such solution is uniformly bounded: Corollary 2. Suppose we are under the hypotheses of Corollary 1.If there exists a solution x * ∈ C([0, +∞[, X) and a constant k > 0 such that for any t ≥ 0 we have |x * (t)| ≤ k, then such solution is unique. .
Fix T 1 = ∆T and observe that, by using Corollary 1, there exists a unique solution in [0, T 1 ].Since x * is a solution of such problem, we have that x * is unique.Suppose we have defined T n−1 such that x * is the unique solution of system (10) Define then T n = T n−1 + ∆T and observe that, since |x * (T n−1 )| ≤ k, by using Corollary 1, there exists a unique solution in [T n−1 , T n ].
By using a change of variables, it is easy to show that . By using such relation, we have that system ( 11) is equivalent to whose unique solution is x(t) = x * (t + T n−1 ) so that x(t) = x * (t) and x * (t) is the unique solution of system (10) in [0, T n ].Since T n → +∞ as n → +∞, we have global uniqueness of limited solutions.

The Fractional M/M/1 Queue
Let us consider again the fractional M/M/1 process N ν (t), t ≥ 0 defined by (3) with state probabilites in (6).
Consider the Hilbert space (l ) be the space of the ν-Hölder continuous functions from [0, T] to l 2 (R).One can rewrite the system (5) in l 2 (R) as follows: where is an infinite tridiagonal matrix with A 0 = (a i,j ) i,j≥0 .Let us show the following: Proof.To show that A 0 is continuous, let us use Schur's test (Theorem 5.2 in [21]).Observe that Moreover, By Schur's test, we have that A 0 is a bounded operator on l 2 and Thus, by Corollary 1, we obtain local existence and uniqueness of the solution of system (5).Global uniqueness can be obtained a posteriori, since the solutions of such system are known.
Let us also observe that the distributions of the inter-arrival times are Mittag-Leffler distributions.To do that, consider the system, for fixed n ≥ 0 which are the state probabilities of a queue with null death rate, fixed birth rate, starting with n customers and with an absorbent state n + 1.Under such assumptions, b ν n+1 (t) is the probability that a customer arrives before t.Moreover, the normalizing condition becomes One can solve the first equation (see Appendix A) to obtain where E ν is the one-parameter Mittag-Leffler function (see Appendix B), and then, by using the normalizing condition, we have In a similar way, let us show that the distributions of the service times are Mittag-Leffler distributions.To show that, consider the system, for fixed n ≥ 0, which are the state probabilities of a queue with null birth rate, fixed death rate, starting with n + 1 customers with an absorbent state n.Under such assumption, d ν n (t) is the probability that a customer is served before t.Moreover, the normalizing condition becomes One can solve the second equation to obtain and then, by using the normalizing condition, we have Moreover, since we know that ∀t ≥ 0 p ν n (t) ≥ 0 and [22]), (p n (t)) n≥0 is uniformly bounded in l 2 (R) and then, by Corollary 2, it is the (global) unique solution of system (5).

Distribution of the Busy Period
We want to determine the probability distribution K ν (t) of the busy period K ν of a fractional M/M/1 queue.To do this, we will follow the lines of the proof given in [1] and [4].Theorem 2. Let K ν be the random variable describing the duration of the busy period of a fractional M/M/1 queue N ν (t) and consider K ν (t) = P(K ν ≤ t).Then, where Proof.Let us first define a queue N ν (t) such that P(N except for the state 0 being an absorbent state.Thus, state probability functions are solution of the following system First, we want to show that, if we consider L ν (t), the inverse of a ν-stable subordinator that is To do that, consider the probability generating From system (15), we know that G ν (z, t) solves the following fractional Cauchy problem: which, for ν = 1, becomes Taking the Laplace transform in Equation ( 17) and using Equation (A1), we have where G ν (z, s) and π ν 0 (s) are Laplace transforms of G ν (z, t) and p ν 0 (t).We know that N ν (t) and then if and only if, by Equation ( 16), Taking the Laplace transform in Equations ( 20) and (21) for n = 0 and by using (see, i.e., Equation (10) in [12]) we know we have to show that and Since Equation ( 17) admits a unique solution, then we only need to show that the right-hand sides of Equations ( 23) and (24) solve Equation (19), that is to say that we have to verify To do that, consider the right-hand side of the previous equation and, recalling that G 1 (z, t) is solution of Equation ( 18 and then, by integrating by parts, we have Equation (25).Now remark that p ν 0 (t) = B ν (t).Thus, we want to determine p ν 0 (t).To do that, let us recall, from [1,4] that from which, explicitly writing I n (2 αβt), we have and then and then, using Equation (26), we have Taking the Laplace transform in Equation ( 27), using Equation ( 22), we have and then integrating Finally, using formula (A2), we have As ν → 1 we obtain, by using that p ν 0 (t) → p 1 0 (t) and then K ν (t) → K 1 (t).

The Fractional M/M/1 Queue with Catastrophes
Let us consider a classical M/M/1 queue with FIFO discipline and subject to catastrophes whose effect is to instantaneously empty the queue [14] and let N 1 ξ (t) be the number of customers in the system at time t with state probabilities Then, the function p 1,ξ n satisfy the following differential-difference equations: where δ n,0 is the Kroeneker delta symbol, D t = d dt , α, β > 0 are the entrance and service rates, respectively, and ξ > 0 is the rate of the catastrophes when the system is not empty.
For ν ∈ (0, 1), we define the fractional M/M/1 queue process with catastrophes as where L ν is an inverse ν-stable subordinator that is independent of N 1 ξ (t), t ≥ 0 (see Section 2).We will show that the state probabilities satisfy the following differential-difference fractional equations: where C 0 D ν t is the Caputo fractional derivative (see Appendix A).In the classical case, catastrophes occur according to a Poisson process with rate ξ if the system is not empty.In our case, consider for a fixed n > 0, which describes the state probabilities of an initially non empty system with null birth and death rate but positive catastrophe rate.In such case, c ν 0 (t) is the probability a catastrophe occurs before time t.Moreover, the normalization property becomes In such case, we can solve the second equation to obtain Using the normalization property, we finally obtain and then the distributions of the inter-occurrence of the catastrophes are Mittag-Leffler distributions.We can conclude that, in the fractional case, catastrophes occur according to a fractional Poisson process ( [10,11,13]) with rate ξ if the system is not empty.Since the operators C 0 D ν t are Caputo fractional derivatives, we expect the stationary behaviour of the queue to be the same as the classic one.Denoting with N 1 ξ the number of customers in the system at the steady state of a classical M/M/1 with catastrophes and defining the state probabilities we can use the results obtained in [15] to observe that where z 1 is the solution of Let us call z 2 the other solution of Equation ( 32) and observe that 0 < z 2 < 1 < z 1 .Some properties coming from such equations that will be useful hereafter are and with i = 1, 2.

Alternative Representation of the Fractional M/M/1 Queue with Catastrophes
We want to obtain an alternative representation of the fractional M/M/1 queue with catastrophes in a way which is similar to Lemma 2.1 in [14].To do that, we firstly need to assure that system (29) admits a unique uniformly bounded solution.To do that, let us write system (29) in the form where p ν,ξ (t) = (p ν,ξ n (t)) n≥0 ∈ C([0, T], l 2 (R)), f (x) = A ξ x + ξ, ξ = (ξ, 0, . . ., 0, . . . ) and is an infinite tridiagonal matrix with A ξ = (a i,j ) i,j≥0 .We need to show the following: Lemma 2. The linear operator A ξ is continuous and A ξ ≤ 2(α + β) + ξ.
Proof.To obtain an estimate of the norm of A ξ , let us use Schur's test.Observe that By Schur's test, we have that A ξ is a bounded operator on l 2 and Observe that, if ξ = 0, the operator A 0 is the same of system (12).Let us also observe that by Corollary 1 locally exists a unique solution.Moreover, if we show that a solution is uniformly bounded, such solution is unique.Now, we are ready to adapt Lemma 2.1 of [14] to the fractional case.Theorem 3. Let N ν (t) be the number of customers in a fractional M/M/1 queue with arrival rate αz 1 and service rate β z 1 such that P( N ν (0) = 0) = 1 and consider N a random variable independent from N ν (t) whose state probabilities q n are defined in Equation (31).Define M ν (t) := min{ N ν (t), N}, t ≥ 0.
Then, the state probabilities of M ν (t) are the unique solutions of (29).
which, by using the definitions of p ν n (t) and q n , becomes Moreover, by using Equation (31), we have and then, substituting Equation (37) in (36), we obtain We want to show that M ν (t) = N ν (t).Since, by definition, p * ,ν n (t) are non-negative and ∑ +∞ n=0 p * ,ν n (t) = 1, they are uniformly bounded in l 2 (R).Thus, we only need to show that p * ,ν (t) = (p * ,ν n (t)) n≥0 solves system (35).The initial conditions are easily verified, so we only need to verify the differential relations.Observe that and then, applying the Caputo derivative operator, we obtain Since p ν 0 (t) is a solution of system ( 5) with rates αz 1 and β z 1 , we have Observe also that so we have After some calculations, we obtain Let us remark that By using Equations ( 32) and (34), we obtain Rewrite now Equation (38) in the form and then apply a Caputo derivative operator to obtain Since p ν n (t) is a solution of system ( 5) with birth rate αz 1 and death rate β z 1 , then we have Remark that, by using Equation (39), Then, recalling that by definition q n−1 = z 1 q n and q n+1 = q n z 1 and doing some calculations, we have Finally, by using Equations (32), ( 33) and (34), we have We have shown that the state probabilities p * ,ν n (t) of M ν (t) are the unique solutions of system (29).Now, we need to show that M ν (t) To do this, consider N 1 (t) a classical M/M/1 queue with arrival rate αz 1 and service rate , N a random variable independent from N ν (t) and N 1 (t) with probability masses q n and finally L ν (t) the inverse of a ν-stable subordinator which is independent from N and N 1 (t).Define also M 1 (t) = min{ N 1 (t), N}.By Lemma 2.1 of [14], we know that M 1 (t) . However, by definition, we know that

State Probabilities for the Fractional M/M/1 with Catastrophes
Since we have defined , where L ν (t) is the inverse of a ν-stable subordinator, which is independent from N 1 ξ (t), we can use such definition and Theorem 3 with the results obtained in [14] to study the state probabilities of N ν ξ (t).In particular, we refer to the formula where By using such formula, we can show the following: Theorem 4. For any t > 0 and n = 0, 1, . . ., we have where C i n,m,r are defined in (41).
and then, by using formula (40), Taking the Laplace transform and using Equation ( 22), we obtain and then, integrating Finally, by using Equation (A2), we obtain  If α > β and ξ → 0 then z 1 = 1 and z 2 = β α .In such case, q n = 0, C 1 n,m,r = 0 and C 2 n,m,r = α n+m β r−n−1 ( m+r m ) r−m m+r .For such case, we have which is not recognizable as a previously obtained formula.This is due to the fact that the formula (which is the one that is obtained from (42) as ν = 1 and α > β, as done in [14]) has no known equivalent in the fractional case.It is also interesting to observe that in [8] another representation of the Laplace transform of p ν n (t) is given in formula 2.40, which is not easily invertible, but has been obtained by using (43) instead of Sharma's representation of p

Distribution of the Busy Period
Let B ν denote the duration of the busy period and B ν (t) = P(B ν ≤ t) be its probability distribution function.Let us observe that, if we pose N ν (0) = 1, then the queue empties within t if and only if a catastrophe occurs within t or otherwise the queue empties without catastrophes within t.Let us remark that, if there is no occurrence of catastrophes, the queue behaves as a fractional M/M/1.Let us define K ν as the duration of a busy period for a fractional M/M/1 queue without catastrophes, Ξ ν the time of first occurrence of a catastrophe for a non empty queue and K ν (t) = P(K ν ≤ t) and Ξ ν (t) = P(Ξ ν ≤ t) their probability distribution functions.Thus, we have Remark 5.If we denote with b ν (t), ξ ν (t) and k ν (t) the probability density functions of B ν , Ξ ν and K ν , we have, by deriving formula (44 which, for ν = 1, is formula (17) of [14].
By using formula (44), we can finally show: Theorem 5. Let B ν be the duration of the busy period of a fractional M/M/1 queue with catastrophes and B ν (t) = P(B ν ≤ t).Then, where C n,m is given in (14).
Proof.Observe that, by formula (30), we have and by formula (13) we also have a closed form of K ν (t).Thus, by using formula (44), we obtain Equation (45).

Distribution of the Time of the First Occurrence of a Catastrophe
We already know that if the queue starts from a non-empty state, then the occurrence of the catastrophes is a Mittag-Leffler distribution.However, we are interested in such distribution as the queue starts being empty.To do that, we will need some auxiliary discrete processes.Theorem 6.Let D ν be the time of first occurrence of a catastrophe as P(N ν (0) = 0) = 1 and let D ν (t) = P(D ν ≤ t).Then, where Proof.Following the lines of [14], let us consider the process N ν (t) with state space {−1, 0, 1, 2, . . .} such that P(N ν (t) = 0) = 1 and posing r n (t) = P(N ν (t) = n), n ≥ −1 as its state probability, we have Let us remark that such process represents our queue until a catastrophe occurs: in such case, instead of emptying the queue, the state of the process becomes −1, which is an absorbent state.With such interpretation, we can easily observe that D ν (t) = r ν −1 (t).
t) be the probability generating function of N ν (t) + 1. Multiplying the third sequence of equations in (47) with z n+1 and then, summing all these equations, we have Now observe that From Theorem 3.1 of [14], we know that and, since we know that N ν (t) , we can use (59) in ( 56) with n = 0 to obtain: Taking the Laplace transform in (60) and using formula (22), we obtain and then integrate (61) Finally, applying the inverse Laplace transform on Equation (61) and using formula (A2), we complete the proof.

Conclusions
Our work focused on the transient behaviour of a fractional M/M/1 queue with catastrophes, deriving formulas for the state probabilities, the distribution of the busy period and the distribution of the time of the first occurrence of a catastrophe.This is a non-Markov generalization of the classical M/M/1 queue with catastrophes, obtained through a time-change.The introduction of fractional dynamics in the equations that master the behaviour of the queue led to a sort of transformation of the time scale.Fractional derivatives are global operators, so the state probabilities preserve memory of their past, eventually slowing down the entire dynamics.Indeed, we can see how Mittag-Leffler functions take place where in the classical case we expected to see exponentials.However, such fractional dynamic seems to affect only the transient behaviour, since we have shown in Remark 2 that the limit behaviour is the same.
The main difficulty that is linked with fractional queues (or in general time-changed queues) is the fact that one has to deal with non-local derivative operators, such as the Caputo derivative, losing Markov property.However, fractional dynamics and fractional processes are gaining attention, due to their wide range of applicability, from physics to finance, from computer science to biology.Moreover, time-changed processes have formed a thriving field of application in mathematical finance. in which the initial condition is given in terms of fractional integrals, while if we use Caputo derivatives we have: C t 0 D ν t x = f (t, x(t)), x(t 0 ) = x 0 , in which the initial condition is related only with the initial value of the function.For such reason, we will prefer adopting Caputo derivatives as fractional derivatives in this paper.
Finally, let us remark that the definition of fractional integral and derivative can be also considered for functions x : [t 0 , t 1 ] ⊆ R → B where B is a Banach space and all the involved integrals are Bochner integrals ( [23]).

Appendix B. Some Special Functions
We recall the definitions of some special functions we use in such text.

Gamma funcion
The Gamma function is defined as: In particular, we have Γ(z + 1) = zΓ(z) and, for z = n ∈ N, Γ(n + 1) = n!.The modified Bessel function ( [24]) of the first kind can be defined by its power series expansion as: One-parameter Mittag-Leffler functions ( [6]) are defined by their power series expansion as: , ν > 0, z ∈ C.
Generalized Mittag-Leffler functions are defined by their power series expansion as: An alternative way to define the Generalized Mittag-Leffler function is: , z ∈ C.
From formula (42), we can easily see that lim t→+∞ p (t) = q n so, as we expected, the steady-state probabilities are the same as the classical ones.For such reason, we can say that the fractional behaviour is influential only in the transient state of the queue.
Remark 2.n In order to determine r ν n (t), we will first show that N ν (t) d = N 1 (L ν (t)) where L ν (t) is the inverse of a ν-stable subordinator which is independent from N 1 .To do that, let us consider N ν (t) + 1 instead of N ν (t).Let us remark that P(N ν (t) + 1