1 Introduction
Matrix Mittag-Leffler (MML) distributions, defined by Albrecher et al. [2, 4], is a wide class of heavy-tailed distributions with some attractive mathematical properties. It is dense in the class of all lifetime distributions, like Phase Type (PH) distribution. This distribution is helpful in simulating a variety of real-world situations, where inter-arrival times of a queueing system have thicker tails. Additionally, the class of MML distributions is a fractional generalisation of a PH distribution (referred to it as a fractional PH distribution [3]). The PH distribution, the Erlang distribution, the exponential distribution, the Mittag-Leffler (ML) distribution, and the fractional Erlang distribution are special cases of this distribution(see, [16, 21, 18, 23, 9, 24, 2, 3]). As introduced by Albrecher et al. [2], PH distribution plays a very important role in MML distribution. A random variable X is said to have a MML distribution if its Laplace transform is given by $\mathbb{E}({e^{-uX}})=\boldsymbol{\pi }{({u^{\alpha }}I-\textbf{T})^{-1}}t$ and is denoted as $X\sim MML(\alpha ,\boldsymbol{\pi },\textbf{T})$, where $(\boldsymbol{\pi },\textbf{T})$ is PH representation and $0\lt \alpha \le 1$. Note that for $\alpha =1$, X follows a PH distribution. PH distributions are dense (weakly converging) on the positive real line, allowing them to accurately resemble a positive distribution. However, their light tail may be a challenge in applications that rely largely on tail behaviour [8]. Fitting heavy-tailed distributions with a PH distribution can result in the need for several stages and the resulting model may still fail to correctly represent the tail behaviour. MML distributions possess advantageous characteristics for representing heavy-tailed phenomena, and they can surpass alternative modelling methods in a noteworthy manner. Moreover, this type of distributions is a specific expansion of the Inhomogeneous Phase Type (IPH) class [2].
Due to the heavy tailed nature of MML distributions, they become very useful in modeling events in many real-life applications in queueing systems, especially when inter-arrival times are heavy-tailed. Heavy-tailed distributions appear naturally in queueing models (see [10, 15, 29, 14]). In the past 100 years, researchers and practitioners have given queueing models a lot of attention, resulting in an extensive collection of papers and books (see, [13, 17, 26, 7, 1, 11, 27]). Recently, researchers have examined fractional queue models that involve catastrophes in the model (see [5, 28, 11, 12, 25]). There is a correlated fractional Erlang queue model given in [6]. M/M/1 queue model employs exponentially distributed inter-arrival and service times for the clients, making it suitable for rigorous mathematical analysis and serving as a useful starting point for various scenarios. As a basic model, it does not include certain characteristics, such as memory preservation, which are frequently needed for analysing more intricate systems. There are multiple instances where memory effects are essential. For this reason, more flexible interarrival-time models are needed in queueing theory. PH distributions and, more generally, Mittag–Leffler-type interarrival times are well suited for this purpose, as they can capture heavy-tailed behavior and dependence structures while still retaining analytical tractability. These features make them realistic and powerful tools for modeling complex arrival processes in modern queueing systems.
In this article, we derive explicit forms of the MML density functions for the truncated queueing systems PH/M/$1\sim n$ and M/PH/$1\sim n$, where n denotes the truncation level at which the queueing process is absorbed. We first study PH/M/1 queues with PH-distributed inter-arrival times, focusing on the exponential and Erlang cases, and obtain explicit expressions for the corresponding MML densities. We then consider M/PH/1 queues with PH-distributed service times and derive analogous formulas when the service-time distribution is exponential or Erlang. These results provide tractable representations of the underlying MML distributions and illustrate how PH-structured arrival and service mechanisms naturally give rise to Mittag-Leffler type behavior in queueing systems. Such density representations may be useful in modeling systems exhibiting heavy-tailed inter-arrival or service-time behavior.
The rest of the paper is organized as follows. In Section 2, we define a PH distribution and a matrix Mittag-Leffler function. We derive the density function of MML random variable for the PH / M /1 $\sim n$ system in Section 3. Section 4 presents the density function of MML random variable for M/PH/1 $\sim n$ system. The article is concluded with some discussions in Section 5.
2 Phase-Type Distributions(PH-distributions) and Matrix Mittag-Leffler function
2.1 Infinitesimal generator
An infinitesimal generator is a square matrix, of order infinite or finite, such that (i) every diagonal element is non-positive,(ii) every off-diagonal element is non-negative,(iii) each row sum is zero.
Infinitesimal generators arise naturally in continuous-time Markov chains (CTMC). For a CTMC with generator Q, the transition probabilities over time $t\ge 0$ are given by the transition semigroup $P(t)={e^{tQ}}$. The $(i,j)$th entry of $P(t)$ represents the probability of being in state j at time t given that the chain started in state i. Thus, the matrix Q characterizes the instantaneous transition rates of the CTMC and determines its entire probabilistic evolution through ${e^{tQ}}$. The diagonal element ${q_{ii}}$ represents the total rate of leaving state i, so the waiting time in state i is exponentially distributed with parameter $-{q_{ii}}$. The off-diagonal element ${q_{ij}}$ for $i\ne j$ gives the rate of transition from state i to state j, and the ratio $\frac{{q_{ij}}}{-{q_{ii}}}$ gives the probability of jumping to state j when a transition occurs.
For example, a 5-order infinitesimal generator, corresponding to a continuous Markov chain with 5-states $\{1,2,3,4,5\}$, is given by(see, [19])
From the definition of a continuous-time Markov chain, the waiting time in each of the states, $\{1,2,3,4,5\}$, is exponentially distributed with parameters $5,\hspace{3.33333pt}10,\hspace{3.33333pt}5,\hspace{3.33333pt}1$ and 0, respectively. State 5 is unique since its waiting time parameter is zero, indicating that there is no state transition after the Markov chain reaches this state. Thereafter, the Markov chain will always be in the state 5, the absorption state.
(1)
\[ Q=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}-5& 0& 1& 0& 4\\ {} 1& -10& 0& 0& 9\\ {} 0& 2& -5& 0& 3\\ {} 0& 0& 1& -1& 0\\ {} 0& 0& 0& 0& 0\end{array}\right).\]Now consider a continuous time Markov chain $I(t)$ having $n+1$ states $\{1,2,3,\dots ,n+1\}$ with its infinitesimal generator
where, T is a PH-generator of order n, that is, T is $n\times n$ invertible matrix such that(i) every diagonal element is negative,(ii) every off-diagonal element is non-negative,(iii) each row sum is non-positive.A PH-generator arises in the context of phase-type (PH) distributions, which model the time until absorption in a finite-state CTMC with one absorbing state. A formal definition of Phase-Type (PH) distributions will be given later in this section.
Since each row sum of matrix Q is zero, we have ${T^{0}}=-T\textbf{e}$, where e is a column vector of order n, having all elements as one. As the total transition rate is zero, one of the states is absorption state, say $n+1$.
Let $I(t)$ denote the continuous-time Markov chain (CTMC) with state space $\{1,2,\dots ,n+1\}$ and infinitesimal generator Q, where the initial distribution is $(\beta ,0)$. Define a random variable
that is, the absorption time of state $n+1$. We say that X is a Phase-Type (PH) random variable with PH representation $(\beta ,T)$ (Neuts [22]).
It is not difficult to check that the distribution of X is given by
For the Markov chain corresponding to Q, given in (1), state 5 is an absorption state. The absorption time of state 5 has a PH-distribution with $m=\hspace{3.33333pt}4$ and
(2)
\[ {F_{X}}(t)=P\{X\le t\}=1-\beta {e^{(Tt)}}\mathbf{e}\equiv 1-\beta \Big({\sum \limits_{k=0}^{\infty }}\frac{{t^{k}}}{k!}{T^{k}}\Big)\mathbf{e},\hspace{3.33333pt}\hspace{3.33333pt}t\ge 0,\]
\[ T=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}-5& 0& 1& 0\\ {} 1& -10& 0& 0\\ {} 0& 2& -5& 0\\ {} 0& 0& 1& -1\end{array}\right),\hspace{3.33333pt}\hspace{3.33333pt}\hspace{3.33333pt}\hspace{3.33333pt}\hspace{3.33333pt}{T^{0}}=\left(\begin{array}{c}4\\ {} 9\\ {} 3\\ {} 0\end{array}\right).\]
2.2 Matrix Mittag–Leffler(MML) function
Let $A\in {\mathbb{C}^{n\times n}}$ (set of $n\times n$ matrices of complex numbers) and, α and β be two complex numbers such that $Re(\alpha )\gt 0$. The MML function has the following equivalent forms:([2, 24])
In this article, we use the Cauchy-integral form of MML function.
-
• Jordan canonical formLet ${\lambda _{1}},{\lambda _{2}},\dots {\lambda _{p}}$ be distinct eigenvalues of A, then A can be expressed in the Jordan canonical form with ${J_{k}}=$ $\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}{\lambda _{k}}& 1& 0& \dots & 0\\ {} 0& {\lambda _{k}}& 1& \dots & 0\\ {} 0& 0& {\lambda _{k}}& \dots & 0\\ {} \vdots & \vdots & \vdots & \dots & \vdots \\ {} 0& 0& 0& \dots & {\lambda _{k}}\end{array}\right)\in {\mathbb{C}^{{m_{k}}\times {m_{k}}}}$,and ${m_{1}}+\cdots +{m_{p}}=n$.The MML function is then expressed as\[ {E_{\alpha ,\beta }}(A)=P\hspace{3.33333pt}\text{diag}({E_{\alpha ,\beta }}({J_{1}}),\dots ,{E_{\alpha ,\beta }}({J_{k}})){P^{-1}},\]where, ${E_{\alpha ,\beta }}({J_{k}})=$ $\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}{E_{\alpha ,\beta }}({\lambda _{k}})& {E_{\alpha ,\beta }^{(1)}}({\lambda _{k}})& \frac{{E_{\alpha ,\beta }^{(2)}}({\lambda _{k}})}{2!}& \dots & \frac{{E_{\alpha ,\beta }^{({m_{k}}-1)}}({\lambda _{k}})}{({m_{k}}-1)!}\\ {} 0& {E_{\alpha ,\beta }}({\lambda _{k}})& {E_{\alpha ,\beta }^{(1)}}({\lambda _{k}})& \dots & \frac{{E_{\alpha ,\beta }^{({m_{k}}-2)}}({\lambda _{k}})}{({m_{k}}-2)!}\\ {} 0& 0& {E_{\alpha ,\beta }}({\lambda _{k}})& \dots & \frac{{E_{\alpha ,\beta }^{({m_{k}}-3)}}({\lambda _{k}})}{({m_{k}}-3)!}\\ {} \vdots & \vdots & \vdots & \dots & \vdots \\ {} 0& 0& 0& \dots & {E_{\alpha ,\beta }}({\lambda _{k}})\end{array}\right),$${E_{\alpha ,\beta }^{(i)}}({\lambda _{k}})$ is ${i^{th}}$ derivative of ${E_{\alpha ,\beta }}({\lambda _{k}}).$
2.3 Matrix Mittag–Leffler distribution
As mentioned earlier, the random variable X is said to have a MML distribution if its Laplace transform is given by
and we write $X\sim MML(\alpha ,\boldsymbol{\pi },\textbf{T})$.
Here, $\boldsymbol{\pi }$ denotes a stochastic row vector and T is a PH-generator.
The density function of X is given by
(see, [2], for more details).
(3)
\[ f(x)={x^{\alpha -1}}\boldsymbol{\pi }{E_{\alpha ,\alpha }}(\textbf{T}{x^{\alpha }})t,\hspace{3.33333pt}\hspace{3.33333pt}x\gt 0.\]We note that the MML distribution may admit a probabilistic representation of the form $X={S_{\alpha }}Y$, where ${S_{\alpha }}$ is a positive α-stable random variable and Y is a PH-distributed random variable as introduced in [2, 4, 3]. Such a representation would facilitate simulation of the MML distribution, since ${S_{\alpha }}$ can be generated using Kanter’s formula and PH distributions can be simulated using standard algorithms. A full derivation of this representation is beyond the scope of this work, but it represents an interesting direction for future research.
2.4 Renewal process
Definition 1 (Renewal Process).
A stochastic process $\{N(t),\hspace{0.1667em}t\ge 0\}$ is called a renewal process if it counts the number of events (renewals) that have occurred by time t, where the interarrival times between consecutive events form an independent and identically distributed (i.i.d.) sequence of non–negative random variables ${\{{X_{n}}\}_{n\ge 1}}$. Let
be the time of the n-th renewal. Then the renewal process is defined by
Renewal processes are widely used to model systems that reset or renew after random time intervals, such as component replacements in reliability theory or recurrent events in stochastic modeling.
3 The PH/M/1 queue
In the PH/M/1 queue customers arrive according to a PH-renewal process with PH- representation $(\beta ,T)$ (see, [19] Section 4.1). For each arriving customer, we initiate a CTMC with generator Q given in Section 2.1, and initial distribution $(\beta ,0)$. The next customer arrival occurs when this CTMC reaches its absorbing state. By concatenating these CTMCs over successive arrivals, we define the phase of the arrival process at time t as ${I_{a}}(t)$. Let $q(t)$ denote the length of the queue at time t.
The following result is due to [19].
Theorem 1.
The stochastic process $\{(q(t),{I_{a}}(t)),t\ge 0\}$ is a continuous time Markov chain with state space $\{0,1,2,...\}\times \{1,2,...,{m_{a}}\}$, and infinitesimal generator
\[ Q=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}T& {T^{0}}\beta \\ {} \mu I& -\mu I+T& {T^{0}}\beta \\ {} & \mu I& -\mu I+T& {T^{0}}\beta \\ {} & & \ddots & \ddots & \ddots \end{array}\right),\]
where, T is PH-generator of order ${m_{a}}$.
In this subsection, we consider a truncated version of the PH/M/1 queue.
3.1 MML distribution of PH/M/1 $\sim n$ system
Here, the notation PH/M/1 $\sim n$ indicates that the underlying PH/M/1 queue is stopped once the queue length reaches level $n+1$ that is, the state $n+1$ is taken to be absorbing. The process evolves like a standard PH/M/1 queue until this threshold is hit. Let $q(t)$ denote the number of customers in the system at time t, and let ${I_{a}}(t)$ be the phase of the arrival process at time t. The pair $(q(t),{I_{a}}(t))$ therefore forms a continuous-time Markov chain on the state space $\{0,1,2,...,n\}\times \{1,2,...,{m_{a}}\}$, with absorption occurring when $q(t)=n+1$.
The PH-generator for $\{(q(t),{I_{a}}(t)),t\ge 0\}$ is defined as
The next result proves the existence of a PH-distributed random variable.
Theorem 2.
Let ${\pi _{1}}={(1,0,0,\dots ,0)_{1\times (n+1)}}$ and, β be a stochastic vector of PH distribution of PH/M/1 $\sim n$ queue. If $\boldsymbol{\pi }=(\beta \textstyle\bigotimes {\pi _{1}})$, where ⨂ is Kronecker product and, T is the PH-generator of PH/M/1 $\sim n$ system. Then, there exist a non-negative random variable X that follows PH-distribution with parameter $(\boldsymbol{\pi },\mathbf{T})$.
In Theorems 3–6, we consider the special case where the PH interarrival distribution is exponential with rate λ. Under this assumption, the PH/M/$1\sim n$ queue reduces to an M/M/$1\sim n$ queue, where interarrival times are exponentially distributed with rate λ and service times are exponentially distributed with rate μ.
Theorem 3.
Let $X\sim MML(\alpha ,\boldsymbol{\pi },\mathbf{T})$, where $(\boldsymbol{\pi },\mathbf{T})$ (with PH-generator T having order 2) is PH-representation of PH/M/1 $\sim 1$. Then, the probability density function of X is given by
\[ f(x)={x^{2\alpha -1}}{\lambda ^{2}}\Big(\frac{{E_{\alpha ,\alpha }}({s_{1}})}{({s_{1}}-{s_{2}})}+\frac{{E_{\alpha ,\alpha }}({s_{2}})}{({s_{2}}-{s_{1}})}\Big)\hspace{3.57777pt}\textit{with}\hspace{3.57777pt}{s_{1}},{s_{2}}=\frac{(-2\lambda -\mu \pm \sqrt{{\mu ^{2}}+4\mu \lambda })}{2}{x^{\alpha }},\]
where, λ is inter-arrival rate and μ is service rate.
Proof.
Note that, for $n=1$,
\[ \textbf{T}=\left(\begin{array}{c@{\hskip10.0pt}c}-\lambda & \lambda \\ {} \mu & -(\mu +\lambda )\end{array}\right).\]
Clearly, $t=-\textbf{Te}={(0,\lambda )^{T}}$.Thus,
\[ {(sI-{x^{\alpha }}\textbf{T})^{-1}}=\frac{1}{{s^{2}}+s(2\lambda +\mu ){x^{\alpha }}+{\lambda ^{2}}{x^{2\alpha }}}\left(\begin{array}{c@{\hskip10.0pt}c}s+\mu {x^{\alpha }}+\lambda {x^{\alpha }}& {x^{\alpha }}\lambda \\ {} {x^{\alpha }}\mu & s+{x^{\alpha }}\lambda \end{array}\right).\]
Using (3), and the Cauchy integral form, we get
\[ f(x)={x^{\alpha -1}}\frac{1}{2\pi i}{\int _{\Gamma }}{E_{\alpha ,\alpha }}(s)\boldsymbol{\pi }{(sI-{x^{\alpha }}\textbf{T})^{-1}}\hspace{3.33333pt}t\hspace{3.33333pt}ds.\]
This implies, using $\boldsymbol{\pi }=(1,0)$,
\[ f(x)={x^{\alpha -1}}\frac{1}{2\pi i}{\int _{\Gamma }}{E_{\alpha ,\alpha }}(s)\frac{({x^{\alpha }}{\lambda ^{2}})}{{s^{2}}+s(2\lambda +\mu ){x^{\alpha }}+{\lambda ^{2}}{x^{2\alpha }}}ds.\]
Using Cauchy’s residue theorem, we get
where
□The following example demonstrates the above result for specific values of parameters.
Example 1.
Consider $\lambda =2$ and $\mu =1$, that is,
and, $t=-\textbf{\textit{Te}}={(0,2)^{T}}$.
In this case, ${s_{1}}=-{x^{\alpha }},\hspace{3.57777pt}\hspace{3.57777pt}{s_{2}}=-4{x^{\alpha }}$
Thus, we get
\[\begin{aligned}{}f(x)& =4{x^{2\alpha -1}}\Big(\frac{{E_{\alpha ,\alpha }}(-{x^{\alpha }})}{3{x^{\alpha }}}+\frac{{E_{\alpha ,\alpha }}(-4{x^{\alpha }})}{-3{x^{\alpha }}}\Big)\\ {} & =\frac{4}{3}{x^{\alpha -1}}\Big({E_{\alpha ,\alpha }}(-{x^{\alpha }})-{E_{\alpha ,\alpha }}(-4{x^{\alpha }})\Big).\end{aligned}\]
To provide numerical support for the explicit MML density derived in Example 1, we compute and plot the density using the closed-form expression given in that example. This numerical illustration confirms the validity and numerical stability of the derived formula and demonstrates its potential applicability in practical queueing scenarios.
The next result is extension of Theorem (3) for underlying PH-generator of order 3.
Theorem 4.
Let $X\sim MML(\alpha ,\boldsymbol{\pi },\mathbf{T})$, where $(\boldsymbol{\pi },\mathbf{T})$ (with PH-generator T having order 3) is PH-representation of PH/M/1 $\sim 2$. Then, the density function of X is given by
\[\begin{array}{c}\displaystyle f(x)=\frac{{x^{\alpha -1}}}{2\pi i}{\int _{\Gamma }}{E_{\alpha ,\alpha }}(s)\frac{({x^{2\alpha }}{\lambda ^{3}})}{\Big({s^{3}}+{s^{2}}(3\lambda +2\mu ){x^{\alpha }}+s(3{\lambda ^{2}}+2\lambda \mu +{\mu ^{2}}){x^{2\alpha }}}ds,\\ {} \displaystyle +{\lambda ^{3}}{x^{3\alpha }}\Big)\end{array}\]
where λ is inter-arrival rate and μ is service rate.
Proof.
Note that, for $n=2$
\[ \textbf{T}=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c}-\lambda & \lambda & 0\\ {} \mu & -(\mu +\lambda )& \lambda \\ {} 0& \mu & -(\mu +\lambda )\end{array}\right).\]
Thus, $t=-\textbf{Te}={(0,0,\lambda )^{T}}$.Using (3), and the Cauchy integral form, we get
\[ f(x)={x^{\alpha -1}}\frac{1}{2\pi i}{\int _{\Gamma }}{E_{\alpha ,\alpha }}(s)\boldsymbol{\pi }{(sI-{x^{\alpha }}\textbf{T})^{-1}}\hspace{3.33333pt}t\hspace{3.33333pt}ds.\]
This implies, using $\boldsymbol{\pi }=(1,0,0)$ and ${(sI-{x^{\alpha }}\textbf{T})^{-1}}$, we get the desired result. □The following corollary is a special case when all the roots of
are distinct.
Corollary 1.
If all the roots of (4) are distinct then the density function has the following form
\[ f(x)={x^{3\alpha -1}}{\lambda ^{3}}{\sum \limits_{k=0}^{2}}\frac{{E_{\alpha ,\alpha }}({s_{k}})}{({s_{k}}-{s_{k+1}})({s_{k}}-{s_{k+2}})},\]
where
\[ {s_{k}}=\frac{-1}{3}\Big((3\lambda +2\mu ){x^{\alpha }}+{\xi ^{k}}C+\frac{({\mu ^{2}}+6\mu \lambda ){x^{2\alpha }}}{{\xi ^{k}}C}\Big),\hspace{3.57777pt}\hspace{3.57777pt}k\in {\mathbb{Z}_{3}},\]
here ${\xi ^{k}}$ is the cube roots of unity and ${\mathbb{Z}_{3}}$ is a group under addition mod 3 and
We consider the following example to demonstrate the above result.
Example 2.
Consider $\lambda =1$ and $\mu =1$, then we have.
\[ \textbf{\textit{T}}=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c}-1& 1& 0\\ {} 1& -2& 1\\ {} 0& 1& -2\end{array}\right).\]
Thus, $t=-\textbf{\textit{Te}}={(0,0,1)^{T}}$.
Using the above corollary, we get
\[\begin{aligned}{}f(x)& ={x^{3\alpha -1}}{\sum \limits_{k=0}^{2}}\frac{{E_{\alpha ,\alpha }}({s_{k}})}{({s_{k}}-{s_{k+1}})({s_{k}}-{s_{k+2}})},\hspace{3.57777pt}\hspace{3.57777pt}\hspace{3.57777pt}k\in {\mathbb{Z}_{3}}\end{aligned}\]
where ${s_{0}}=-3.25{x^{\alpha }},\hspace{3.57777pt}{s_{1}}=-1.55{x^{\alpha }},\hspace{3.57777pt}{s_{2}}=-0.2{x^{\alpha }}$. Therefore,
The next result is extension of Theorem (4) when the underlying PH-generator has order 4.
Theorem 5.
Let $X\sim MML(\alpha ,\boldsymbol{\pi },\mathbf{T})$, where $(\boldsymbol{\pi },\mathbf{T})$ (with PH-generator T having order 4) is PH-representation of PH/M/1 $\sim 3$. Then, the density function of X is given by
\[\begin{array}{c}\displaystyle f(x)=\frac{{x^{\alpha -1}}}{2\pi i}{\int _{\Gamma }}\frac{{E_{\alpha ,\alpha }}(s)\hspace{3.57777pt}\hspace{3.57777pt}({x^{3\alpha }}{\lambda ^{4}})}{\Big({s^{4}}+{s^{3}}(4\lambda +3\mu ){x^{\alpha }}+{s^{2}}(6{\lambda ^{2}}+6\lambda \mu +3{\mu ^{2}}){x^{2\alpha }}}ds,\\ {} \displaystyle +s(4{\lambda ^{3}}+3{\lambda ^{2}}\mu +2\lambda {\mu ^{2}}+{\mu ^{3}}){x^{3\alpha }}+{\lambda ^{4}}{x^{4\alpha }}\Big)\end{array}\]
where λ is arrival rate and μ is service rate.
Proof.
Note that, for $n=3$,
\[ \textbf{T}=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}-\lambda & \lambda & 0& 0\\ {} \mu & -(\mu +\lambda )& \lambda & 0\\ {} 0& \mu & -(\mu +\lambda )& \lambda \\ {} 0& 0& \mu & -(\mu +\lambda )\end{array}\right).\]
Thus, $t=-\textbf{Te}={(0,0,0,\lambda )^{T}}$.Using (3), and the Cauchy integral form, we have
\[ f(x)={x^{\alpha -1}}\frac{1}{2\pi i}{\int _{\Gamma }}{E_{\alpha ,\alpha }}(s)\boldsymbol{\pi }{(sI-{x^{\alpha }}\textbf{T})^{-1}}\hspace{3.33333pt}t\hspace{3.33333pt}ds.\]
Using $\boldsymbol{\pi }=(1,0,0,0)$ and ${(sI-{x^{\alpha }}\textbf{T})^{-1}}$, we get the desired result. □The following corollary is a special case when all the roots of
are distinct.
Corollary 2.
If all the roots of (5) are distinct then the density function has the following form
\[ f(x)={x^{4\alpha -1}}{\lambda ^{4}}{\sum \limits_{k=0}^{3}}\frac{{E_{\alpha ,\alpha }}({s_{k}})}{({s_{k}}-{s_{k+1}})({s_{k}}-{s_{k+2}})({s_{k}}-{s_{k+3}})}\hspace{3.57777pt},\hspace{3.57777pt}\hspace{3.57777pt}\hspace{3.57777pt}k\in {\mathbb{Z}_{4}},\]
where ${\mathbb{Z}_{4}}$ is a group under addition mod 4 and
\[ {s_{0}},{s_{1}}=-\Big(\frac{4\lambda +3\mu }{4}\Big){x^{\alpha }}-S\pm \frac{1}{2}\sqrt{-4{S^{2}}-2p+\frac{q}{S}}\]
\[ {s_{2}},{s_{3}}=-\Big(\frac{4\lambda +3\mu }{4}\Big){x^{\alpha }}+S\pm \frac{1}{2}\sqrt{-4{S^{2}}-2p-\frac{q}{S}}.\]
Here,
\[ p=\frac{-3\mu }{8}(8\lambda +\mu ){x^{2\alpha }},\hspace{3.57777pt}\hspace{3.57777pt}\hspace{3.57777pt}q=\frac{{\mu ^{2}}}{8}(4\lambda -\mu ){x^{3\alpha }},\]
and
Example 3.
Consider $\lambda =1$ and $\mu =1$, then we have.
\[ \textbf{\textit{T}}=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}-1& 1& 0& 0\\ {} 1& -2& 1& 0\\ {} 0& 1& -2& 1\\ {} 0& 0& 1& -2\end{array}\right).\]
Thus, $t=-\textbf{\textit{Te}}={(0,0,0,1)^{T}}$.
using the above corollary, we get
\[\begin{aligned}{}f(x)& ={x^{4\alpha -1}}{\sum \limits_{k=0}^{3}}\frac{{E_{\alpha ,\alpha }}({s_{k}})}{({s_{k}}-{s_{k+1}})({s_{k}}-{s_{k+2}})({s_{k}}-{s_{k+3}})}\hspace{3.57777pt},\hspace{3.57777pt}\hspace{3.57777pt}\hspace{3.57777pt}k\in {\mathbb{Z}_{4}}\end{aligned}\]
where,
Then
In Theorem (5), we assume that the PH-generator has order 4, the same can be generalized for PH- generator having order n. For $n\ge 5$, obtaining a closed-form expression for the density is generally impossible, since the computation requires solving a polynomial of degree n, and polynomials of degree $n\ge 5$ do not admit closed-form roots in radicals (Abel–Ruffini theorem). Nevertheless, studying the polynomial appearing in the denominator of the Cauchy integral representation of the MML function may still reveal useful structural properties and represents an interesting direction for further investigation. We obtain the density for MML-random variable corresponding to n-order PH-generator when $\mu =0$.
Theorem 6.
Consider PH/M/1 $\sim (n-1)$ with $\mu =0$. Then, the following random variables follow same distribution with common density function(a) $X\sim MML(\alpha ,{\pi _{1}},{T_{1}})$, and (b) $X\sim MML(\alpha ,{\pi _{2}},{T_{2}})$,where $({\pi _{1}},{T_{1}})$ is the PH-representation, given in Theorem (2), of exponential interarrival time with rate λ and $({\pi _{2}},{T_{2}})$ is the PH-representation of Erlang distribution with shape parameter n and rate λ.
Proof.
Note that
\[ {T_{1}}={T_{2}}=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}-\lambda & \lambda & 0& \dots & 0& 0\\ {} 0& -\lambda & \lambda & \dots & 0& 0\\ {} 0& 0& -\lambda & \dots & 0& 0\\ {} \vdots & \vdots & \vdots & \dots & \vdots & \vdots \\ {} 0& 0& 0& \dots & -\lambda & \lambda \\ {} 0& 0& 0& \dots & 0& -\lambda \end{array}\right)\]
(see, [2] and Example 3.1 for ${T_{2}}$)So far, we used PH-representation of exponential distribution. We now replace the assumption of exponential by Erlang and obtain the density of MML random variable.
Theorem 7.
Let $X\sim MML(\alpha ,\boldsymbol{\pi },\mathbf{T})$, where $(\boldsymbol{\pi },\mathbf{T})$ is PH-representation of PH/M/1 $\sim 1$, where the interarrival times follow an Erlang-2 distribution with rate λ. Then, the density function of X is given by
(6)
\[\begin{array}{cc}& \displaystyle f(x)=\frac{{x^{\alpha -1}}}{2\pi i}{\int _{\Gamma }}{E_{\alpha ,\alpha }}(s)\frac{({x^{3\alpha }}{\lambda ^{4}})\hspace{3.57777pt}\hspace{3.57777pt}ds}{\Big({s^{4}}+{s^{3}}(4\lambda +2\mu ){x^{\alpha }}+{s^{2}}(6{\lambda ^{2}}+6\lambda \mu +{\mu ^{2}}){x^{2\alpha }}}ds.\\ {} & \displaystyle +s(4{\lambda ^{3}}+4{\lambda ^{2}}\mu +2\lambda {\mu ^{2}}){x^{3\alpha }}+{\lambda ^{4}}{x^{4\alpha }}\Big)\end{array}\]Proof.
Note that
\[ \textbf{T}=\left(\begin{array}{c@{\hskip10.0pt}c}T& {T^{0}}\beta \\ {} \mu I& -\mu I+T\end{array}\right).\]
Using the PH-representation of Erlang-2 distribution, we have
\[ \textbf{T}=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}-\lambda & \lambda & 0& 0\\ {} 0& -\lambda & \lambda & 0\\ {} \mu & 0& -(\mu +\lambda )& \lambda \\ {} 0& \mu & 0& -(\mu +\lambda )\end{array}\right).\]
Thus, $t=-\textbf{Te}={(0,0,0,\lambda )^{T}}$.Using (3), and the Cauchy integral form, we get
\[ f(x)={x^{\alpha -1}}\frac{1}{2\pi i}{\int _{\Gamma }}{E_{\alpha ,\alpha }}(s)\boldsymbol{\pi }{(sI-{x^{\alpha }}\textbf{T})^{-1}}\hspace{3.33333pt}t\hspace{3.33333pt}ds.\]
The result follows by using $\boldsymbol{\pi }=(1,0,0,0)$ and ${(sI-{x^{\alpha }}\textbf{T})^{-1}}$. □The following corollary gives an explicit form of density, when the roots of
are distinct.
Corollary 3.
If all the roots of (7) are distinct then the density function (6) has the form
\[ f(x)={x^{4\alpha -1}}{\lambda ^{4}}{\sum \limits_{k=0}^{3}}\frac{{E_{\alpha ,\alpha }}({s_{k}})}{({s_{k}}-{s_{k+1}})({s_{k}}-{s_{k+2}})({s_{k}}-{s_{k+3}})}\hspace{3.57777pt},\hspace{3.57777pt}\hspace{3.57777pt}\hspace{3.57777pt}k\in {\mathbb{Z}_{4}},\]
where, ${\mathbb{Z}_{4}}$ is a group under addition mod 4 and
\[ {s_{0}},{s_{1}}=-\Big(\frac{4\lambda +2\mu }{4}\Big){x^{\alpha }}-S\pm \frac{1}{2}\sqrt{-4{S^{2}}-2p+\frac{q}{S}},\]
\[ {s_{2}},{s_{3}}=-\Big(\frac{4\lambda +2\mu }{4}\Big){x^{\alpha }}+S\pm \frac{1}{2}\sqrt{-4{S^{2}}-2p-\frac{q}{S}},\]
where $p=\frac{-{\mu ^{2}}}{2}{x^{2\alpha }},\hspace{3.57777pt}\hspace{3.57777pt}\hspace{3.57777pt}q=-2\lambda {\mu ^{2}}{x^{3\alpha }},\hspace{3.57777pt}\hspace{3.57777pt}S=\frac{1}{2}\sqrt{\frac{-2}{3}p+\frac{1}{3}(Q+\frac{{\mu ^{4}}{x^{4\alpha }}}{Q})}$ and $Q={x^{2\alpha }}\sqrt[3]{54{\mu ^{2}}{\lambda ^{4}}+{\mu ^{6}}+\sqrt{2916{\lambda ^{8}}{\mu ^{4}}+108{\mu ^{8}}{\lambda ^{4}}}}$
Example 4.
Consider $\lambda =\mu =1$, then we have
\[ \textbf{\textit{T}}=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}-1& 1& 0& 0\\ {} 0& -1& 1& 0\\ {} 1& 0& -2& 1\\ {} 0& 1& 0& -2\end{array}\right).\]
This implies $t=-\textbf{\textit{Te}}={(0,0,0,1)^{T}}$.
Using above corollary and $\boldsymbol{\pi }=(1,0,0,0)$, we get
\[ f(x)={x^{4\alpha -1}}{\sum \limits_{k=0}^{3}}\frac{{E_{\alpha ,\alpha }}({s_{k}})}{({s_{k}}-{s_{k+1}})({s_{k}}-{s_{k+2}})({s_{k}}-{s_{k+3}})}\hspace{3.57777pt},\hspace{3.57777pt}\hspace{3.57777pt}\hspace{3.57777pt}k\in {\mathbb{Z}_{4}}\]
where
\[ {s_{0}}=-0.12{x^{\alpha }},\hspace{3.57777pt}\hspace{3.57777pt}{s_{1}}=-1.47{x^{\alpha }},\hspace{3.57777pt}\hspace{3.57777pt}{s_{2}},{s_{3}}=(-2.21\hspace{3.57777pt}\pm \hspace{3.57777pt}i\hspace{3.57777pt}0.98){x^{\alpha }}.\]
Then
\[\begin{array}{c}\displaystyle f(x)={x^{\alpha -1}}\Big(\frac{{E_{\alpha ,\alpha }}(-0.12{x^{\alpha }})}{7.19}-\frac{{E_{\alpha ,\alpha }}(-1.47{x^{\alpha }})}{2.04}+\frac{{E_{\alpha ,\alpha }}((-2.21+.98i){x^{\alpha }})}{5.44+1.15i}\\ {} \displaystyle +\frac{{E_{\alpha ,\alpha }}((-2.21-.98i){x^{\alpha }})}{5.44-1.15i}\Big).\end{array}\]
In Theorem (7), we assume that the number of persons in the queue is 1 and inter-arrival time follows Erlang-2 distribution. The same can be generalized for n persons in the queue and Erlang-n inter-arrival time. However, it is difficult to obtain a closed form expression for the density.
Next we will derive the explicit form of MML distribution for queue in which service time follows a PH distribution.
4 The M/PH/1 queue
We assume that customers arrive with arrival process being Poisson with parameter λ and service time has a common PH-distribution with service rate μ, and PH-representation $(\gamma ,S)$, with S as a PH-generator of order ${m_{s}}$.
Let ${({A_{n}})_{n\ge 1}}$ be the sequence of arrival times. We define the sequence of departure times ${({D_{n}})_{n\ge 1}}$ together with the service phase process ${I_{s}}(t)$ as follows. We set ${I_{s}}(t)=0$ for $t\lt {A_{1}}$. For $t\ge {A_{1}}$, the process ${I_{s}}(t)$ evolves as a continuous-time Markov chain with generator S and initial distribution $(\gamma ,0)$, and we define ${D_{1}}$ as its time of absorption. Suppose that ${D_{n}}$ has been defined. If ${D_{n}}\lt {A_{n+1}}$, then we set ${I_{s}}(t)=0$ for $t\in [{D_{n}},{A_{n+1}})$. At time $t={A_{n+1}}$, the process ${I_{s}}(t)$ restarts as a continuous-time Markov chain with generator S and initial distribution $(\gamma ,0)$, and ${D_{n+1}}$ is defined as its time of absorption. If instead ${D_{n}}\ge {A_{n+1}}$, then ${I_{s}}(t)$ continues evolving as the same Markov chain for $t\ge {D_{n}}$, and ${D_{n+1}}$ is again defined as its time of absorption. In this way, we obtain a process ${I_{s}}(t)$, $t\ge 0$, which is equal to 0 when the system is empty and describes the service phase otherwise. Moreover, the service times
are PH-distributed with PH representation $(\gamma ,S)$.
Consider the stochastic process $\{(q(t),{I_{s}}(t)),t\ge 0\}$. Because the arrival process is a Poisson process, the inter-arrival time t has exponential distribution with rate λ. Furthermore, because the service phase ${I_{s}}(t)$ is recorded, it is possible to calculate the probability of service time completion. Then, it is easy to observe that $\{(q(t),{I_{s}}(t)),t\ge 0\}$ is a continuous time Markov chain with state space $\{0,1,2,...\}\times \{1,2,...,{m_{s}}\}\cup \{0\}$ having infinitesimal generator
\[ Q=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}-\lambda & \lambda \gamma \\ {} {S^{0}}& -\lambda I+S& \lambda I\\ {} & {S^{0}}\gamma & -\lambda I+S& \lambda I\\ {} & & \ddots & \ddots & \ddots \end{array}\right),\]
where, ${S^{0}}=-S\textbf{e}$. The process $\{(q(t),{I_{s}}(t)),t\ge 0\}$ is a QBD (Quasi Birth and Death) process, (see, [19, 20]).4.1 MML distribution for M/PH/1 $\sim n$ system
The notation M/PH/1 $\sim n$ is defined analogously to the PH/M/1 $\sim n$ case in Section 3.1 with state $n+1$ taken as absorbing.
Note that the PH-generator, for $\{(q(t),{I_{s}}(t)),t\ge 0\}$, is defined as
The next result proves the existence of a PH-distributed random variable.
Theorem 8.
Consider M/PH/1 $\sim n$ system with row vector $\boldsymbol{\pi }=(\gamma ,0,0,\dots ,0)$ and PH-generator T. There exists a random variable X which follows PH-distribution with representation $(\boldsymbol{\pi },\mathbf{T})$.
We use Theorem (8) and obtain the distribution of a MML random variable for Erlang underlying PH-distribution. We first obtain the density function of MML random variable, when there is one person in the queue.
Theorem 9.
Let $X\sim MML(\alpha ,\boldsymbol{\pi },\mathbf{T})$, where $(\boldsymbol{\pi },\mathbf{T})$ is PH-representation of M/PH/1 $\sim 1$ system. If the service time follows PH-distribution with PH-representation of an Erlang-2 distribution having service rate μ, that is, $\gamma =(1,0)$. Then, the density function of X is given by
\[\begin{array}{c}\displaystyle f(x)=\frac{{x^{\alpha -1}}}{2\pi i}{\int _{\Gamma }}{E_{\alpha ,\alpha }}(s)\frac{s\lambda {x^{\alpha }}+({\lambda ^{2}}+2\mu \lambda ){x^{2\alpha }}}{\Big({s^{3}}+{s^{2}}(3\lambda +2\mu ){x^{\alpha }}+s(3{\lambda ^{2}}+{\mu ^{2}}+4\mu \lambda ){x^{2\alpha }}}ds,\\ {} \displaystyle +({\lambda ^{3}}+2{\lambda ^{2}}\mu ){x^{3\alpha }}\Big)\end{array}\]
where λ is the arrival rate.
Proof.
Note that
Thus, $t=-\textbf{Te}={(0,\lambda ,\lambda )^{T}}$. We have
\[ \textbf{T}=\left(\begin{array}{c@{\hskip10.0pt}c}-\lambda & \lambda \gamma \\ {} {S^{0}}& -\lambda I+S\end{array}\right),\]
or
(8)
\[ \textbf{T}=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c}-\lambda & \lambda & 0\\ {} 0& -\mu -\lambda & \mu \\ {} \mu & 0& -\mu -\lambda \end{array}\right).\]
\[ {(sI-{x^{\alpha }}\textbf{T})^{-1}}=\frac{1}{\Big({s^{3}}+{s^{2}}(3\lambda +2\mu ){x^{\alpha }}+s(3{\lambda ^{2}}+{\mu ^{2}}+4\mu \lambda ){x^{2\alpha }}+({\lambda ^{3}}+2{\lambda ^{2}}\mu ){x^{3\alpha }}\Big)}\]
\[ \times \left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c}A& B& C\\ {} D& E& F\\ {} G& H& I\end{array}\right)\]
where $A={s^{2}}+s(2\mu +2\lambda ){x^{\alpha }}+({\lambda ^{2}}+{\mu ^{2}}+2\mu \lambda ){x^{2\alpha }},\hspace{3.33333pt}B=s\lambda {x^{\alpha }}+({\lambda ^{2}}+\mu \lambda ){x^{2\alpha }},\hspace{3.33333pt}C=\mu \lambda {x^{2\alpha }},\hspace{3.33333pt}D={\mu ^{2}}{x^{2\alpha }},\hspace{3.33333pt}E={s^{2}}+s(2\lambda +\mu ){x^{\alpha }}+({\lambda ^{2}}+\mu \lambda ){x^{2\alpha }},\hspace{3.33333pt}F=s\mu {x^{\alpha }}+\mu \lambda {x^{2\alpha }},\hspace{3.33333pt}G=s\mu {x^{\alpha }}+({\mu ^{2}}+\mu \lambda ){x^{2\alpha }},\hspace{3.33333pt}H=\mu \lambda {x^{2\alpha }},\hspace{3.33333pt}\text{and},\hspace{3.33333pt}I={s^{2}}+s(2\lambda +\mu ){x^{\alpha }}+({\lambda ^{2}}+\mu \lambda ){x^{2\alpha }}$
Using (3), and the Cauchy integral form, we have
\[ f(x)={x^{\alpha -1}}\frac{1}{2\pi i}{\int _{\Gamma }}{E_{\alpha ,\alpha }}(s)\boldsymbol{\pi }{(sI-{x^{\alpha }}\textbf{T})^{-1}}\hspace{3.33333pt}t\hspace{3.33333pt}ds.\]
The result follows by substituting ${(sI-{x^{\alpha }}\textbf{T})^{-1}}$ and $\boldsymbol{\pi }=(1,0,0)$ in the above integration. □The following corollary is a special case of the above theorem when all the roots of
are distinct.
Corollary 4.
If all the roots of (9) are distinct then the density function of X is given by
\[ f(x)={x^{2\alpha -1}}\lambda {\sum \limits_{k=0}^{2}}\frac{({s_{k}}+(\lambda +2\mu ){x^{\alpha }}){E_{\alpha ,\alpha }}({s_{k}})}{({s_{k}}-{s_{k+1}})({s_{k}}-{s_{k+2}})},\]
where
\[ {s_{k}}=\frac{-1}{3}\Big((3\lambda +2\mu ){x^{\alpha }}+{\xi ^{k}}C+\frac{{\mu ^{2}}{x^{2\alpha }}}{{\xi ^{k}}C}\Big),\hspace{3.57777pt}\hspace{3.57777pt}k\in {\mathbb{Z}_{3}}.\]
Here, ${\xi ^{k}}$ is the cube roots of unity and ${\mathbb{Z}_{3}}$ is a group under addition mod 3, and,
The following example demonstrate the above result for specific values of parameters.
Example 5.
Consider $\lambda =2$ and $\mu =2$. From (8), we have
\[ \textbf{\textit{T}}=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c}-2& 2& 0\\ {} 0& -4& 2\\ {} 2& 0& -4\end{array}\right).\]
Thus, $t=-\textbf{\textit{Te}}={(0,2,2)^{T}}$.
Using above corollary and $\boldsymbol{\pi }=(1,0,0)$ we get
\[\begin{aligned}{}f(x)& ={x^{\alpha -1}}\frac{1}{2\pi i}{\int _{\gamma }}{E_{\alpha ,\alpha }}(s)\frac{2{x^{\alpha }}s+12{x^{2\alpha }}}{{s^{3}}+10{x^{\alpha }}{s^{2}}+32{x^{2\alpha }}s+32{x^{3\alpha }}}ds\\ {} & =2{x^{\alpha -1}}\Big({E_{\alpha ,\alpha }}(-2{x^{\alpha }})-{E_{\alpha ,\alpha }}(-4{x^{\alpha }})-{E^{\prime }_{\alpha ,\alpha }}(-4{x^{\alpha }})\Big),\end{aligned}\]
where ${E^{\prime }_{\alpha ,\alpha }}(-4{x^{\alpha }})$ is the derivative of ${E_{\alpha ,\alpha }}(-4{x^{\alpha }})$.
The next result is extension of Theorem (9), where the underlying PH-distribution is Erlang-3 distribution.
Theorem 10.
Let $X\sim MML(\alpha ,\boldsymbol{\pi },\mathbf{T})$, where $(\boldsymbol{\pi },\mathbf{T})$ is PH-representation of M/PH/1 $\sim 1$ system. If the service time follows PH-distribution with PH-representation of an Erlang-3 distribution with service rate μ, that is, $\gamma =(1,0,0)$. Then, the density function of X is given by
\[ f(x)=\frac{{x^{\alpha -1}}}{2\pi i}\hspace{-0.1667em}\hspace{-0.1667em}{\int _{\Gamma }}\frac{{E_{\alpha ,\alpha }}(s)\Big(\lambda {x^{\alpha }}{s^{2}}+s(2{\lambda ^{2}}+3\lambda \mu ){x^{2\alpha }}+({\lambda ^{3}}+3{\lambda ^{2}}\mu +3\lambda {\mu ^{2}}){x^{3\alpha }}\Big)ds}{{\textstyle\textstyle\sum _{j=0}^{2}}{\textstyle\textstyle\sum _{i=0}^{4-j}}\left(\genfrac{}{}{0.0pt}{}{3}{j}\right)\left(\genfrac{}{}{0.0pt}{}{4-j}{i}\right){\lambda ^{i}}{\mu ^{j}}{x^{(i+j)\alpha }}{s^{4-j-i}}\hspace{2.83862pt}+\hspace{2.83862pt}{\mu ^{3}}{x^{3\alpha }}s,}\]
where λ is the arrival rate.
Proof.
Note that we have
Then $t=-\textbf{Te}={(0,\lambda ,\lambda ,\lambda )^{T}}.$
\[ \textbf{T}=\left(\begin{array}{c@{\hskip10.0pt}c}-\lambda & \lambda \gamma \\ {} {S^{0}}& -\lambda I+S\end{array}\right),\]
that is,
(10)
\[ \textbf{T}=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}-\lambda & \lambda & 0& 0\\ {} 0& -(\mu +\lambda )& \mu & 0\\ {} 0& 0& -(\mu +\lambda )& \mu \\ {} \mu & 0& 0& -(\mu +\lambda )\end{array}\right).\]Using (3), and the Cauchy integral form, we have
\[ f(x)={x^{\alpha -1}}\frac{1}{2\pi i}{\int _{\Gamma }}{E_{\alpha ,\alpha }}(s)\boldsymbol{\pi }{(sI-{x^{\alpha }}\textbf{T})^{-1}}\hspace{3.33333pt}t\hspace{3.33333pt}ds.\]
The result follows by using $\boldsymbol{\pi }=(1,0,0,0)$ and substituting the value of ${(sI-{x^{\alpha }}\textbf{T})^{-1}}$ for T given in (10). □The following corollary is a special case when all the roots of
are distinct.
Corollary 5.
If all the roots of (11) are distinct then the density function has the form
\[ f(x)={x^{2\alpha -1}}\lambda {\sum \limits_{k=0}^{3}}\hspace{-0.1667em}\frac{\Big({s_{k}^{2}}+{s_{k}}(2\lambda +3\mu ){x^{\alpha }}+({\lambda ^{2}}+3\lambda \mu +3{\mu ^{2}}){x^{2\alpha }}\Big){E_{\alpha ,\alpha }}({s_{k}})}{({s_{k}}-{s_{k+1}})({s_{k}}-{s_{k+2}})({s_{k}}-{s_{k+3}})}\hspace{3.57777pt},\hspace{3.57777pt}\hspace{3.57777pt}k\hspace{-0.1667em}\in {\mathbb{Z}_{4}},\]
where ${\mathbb{Z}_{4}}$ is a group under addition mod 4, and
\[ {s_{0}},{s_{1}}=-\Big(\frac{4\lambda +3\mu }{4}\Big){x^{\alpha }}-S\pm \frac{1}{2}\sqrt{-4{S^{2}}-2p+\frac{q}{S}}\]
\[ {s_{2}},{s_{3}}=-\Big(\frac{4\lambda +3\mu }{4}\Big){x^{\alpha }}+S\pm \frac{1}{2}\sqrt{-4{S^{2}}-2p-\frac{q}{S}},\]
where
and
The following example illustrates the above result for specific values of parameters.
Example 6.
Consider $\lambda =1$ and $\mu =1$. We have, from (10),
\[ \textbf{\textit{T}}=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}-1& 1& 0& 0\\ {} 0& -2& 1& 0\\ {} 0& 0& -2& 1\\ {} 1& 0& 0& -2\end{array}\right).\]
Thus $t=-\textbf{\textit{Te}}={(0,1,1,1)^{T}}$.
Now
\[ f(x)={x^{2\alpha -1}}{\sum \limits_{k=0}^{3}}\frac{({s_{k}^{2}}+5{x^{\alpha }}{s_{k}}+7{x^{2\alpha }}){E_{\alpha ,\alpha }}({s_{k}})}{({s_{k}}-{s_{k+1}})({s_{k}}-{s_{k+2}})({s_{k}}-{s_{k+3}})}\hspace{3.57777pt},\hspace{3.57777pt}\hspace{3.57777pt}\hspace{3.57777pt}k\in {\mathbb{Z}_{4}},\]
Where
\[ {s_{0}}-0.62{x^{\alpha }},\hspace{3.57777pt}\hspace{3.57777pt}{s_{1}}=-2.82{x^{\alpha }},\hspace{3.57777pt}\hspace{3.57777pt}{s_{2}},{s_{3}}=(-1.78\pm 0.91\hspace{3.57777pt}i){x^{\alpha }}.\]
Then
\[\begin{array}{c}\displaystyle f(x)={x^{\alpha -1}}\Big(\frac{4.28\hspace{3.57777pt}{E_{\alpha ,\alpha }}(-0.62{x^{\alpha }})}{4.78126}-\frac{0.85\hspace{3.57777pt}{E_{\alpha ,\alpha }}(-2.82{x^{\alpha }})}{4.2416}\\ {} \displaystyle +\frac{(0.44+1.31i){E_{\alpha ,\alpha }}((-1.78+.91i){x^{\alpha }})}{0.2+3.7i}+\frac{(0.44-1.31i){E_{\alpha ,\alpha }}((-1.78-.91i){x^{\alpha }})}{0.2-3.7i}\Big)\end{array}\]
The PH-distribution used in the above Theorem (10) is Erlang-3, we next generalize the result for Erlang-n distribution.
Theorem 11.
Let $X\sim MML(\alpha ,\boldsymbol{\pi },\mathbf{T})$, where $(\boldsymbol{\pi },\mathbf{T})$ is PH-representation of M/PH/1 $\sim 1$ system. If the service time follows PH-distribution with PH-representation of an Erlang-n distribution with service rate μ, that is, $\gamma =(1,0,0,\dots ,0)$. Then the density function of X is given by
\[ f(x)=\frac{{x^{\alpha -1}}}{2\pi i}\hspace{-0.1667em}\hspace{-0.1667em}{\int _{\Gamma }}{E_{\alpha ,\alpha }}(s)\frac{{\textstyle\textstyle\sum _{j=0}^{n-1}}{\textstyle\textstyle\sum _{i=0}^{n-(j+1)}}\left(\genfrac{}{}{0.0pt}{}{n}{j}\right)\left(\genfrac{}{}{0.0pt}{}{n-(j+1)}{i}\right){\lambda ^{i+1}}{\mu ^{j}}{x^{(i+j+1)\alpha }}{s^{n-(j+1)-i}}ds}{{\textstyle\textstyle\sum _{j=0}^{n-1}}{\textstyle\textstyle\sum _{i=0}^{n+1-j}}\left(\genfrac{}{}{0.0pt}{}{n}{j}\right)\left(\genfrac{}{}{0.0pt}{}{n+1-j}{i}\right){\lambda ^{i}}{\mu ^{j}}{x^{(i+j)\alpha }}{s^{n+1-j-i}}+\left(\genfrac{}{}{0.0pt}{}{n}{n}\right){\mu ^{n}}{x^{n\alpha }}s},\]
where λ is the arrival rate.
Proof.
Note that,
\[ \textbf{T}=\left(\begin{array}{c@{\hskip10.0pt}c}-\lambda & \lambda \gamma \\ {} {S^{0}}& -\lambda I+S\end{array}\right),\]
that is,
\[ \mathbf{T}=\left(\begin{array}{c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c@{\hskip10.0pt}c}-\lambda & \lambda & 0& \dots & 0& 0\\ {} 0& -\lambda -\mu & \mu & \dots & 0& 0\\ {} 0& 0& -\lambda -\mu & \dots & 0& 0\\ {} \vdots & \vdots & \vdots & \dots & \vdots & \vdots \\ {} 0& 0& 0& \dots & -\lambda -\mu & \mu \\ {} \mu & 0& 0& \dots & 0& -\lambda -\mu \end{array}\right).\]
Thus, $t=-\textbf{Te}={(0,\lambda ,\lambda ,\dots ,\lambda )_{n+1\times 1}^{T}}.$
Using (3), and the Cauchy integral form, we have
\[ f(x)={x^{\alpha -1}}\frac{1}{2\pi i}{\int _{\Gamma }}{E_{\alpha ,\alpha }}(s)\boldsymbol{\pi }{(sI-{x^{\alpha }}\textbf{T})^{-1}}\hspace{3.33333pt}t\hspace{3.33333pt}ds.\]
The result now follows by using $\boldsymbol{\pi }={(1,0,0,\dots ,0)_{1\times n+1}}$ and substituting the value of ${(sI-{x^{\alpha }}\textbf{T})^{-1}}$ in the above integral. □Following corollary is a special case of the above theorem when all the roots of
are distinct.
(12)
\[ {\sum \limits_{j=0}^{n-1}}{\sum \limits_{i=0}^{n+1-j}}\left(\genfrac{}{}{0.0pt}{}{n}{j}\right)\left(\genfrac{}{}{0.0pt}{}{n+1-j}{i}\right){\lambda ^{i}}{\mu ^{j}}{x^{(i+j)\alpha }}{s^{n+1-j-i}}+\left(\genfrac{}{}{0.0pt}{}{n}{n}\right){\mu ^{n}}{x^{n\alpha }}s=0\]Corollary 6.
If all the roots of (12) are distinct, then the density function is given by
\[ f(x)={x^{2\alpha -1}}\lambda {\sum \limits_{k=0}^{n}}\frac{\Big({\textstyle\textstyle\sum _{j=0}^{n-1}}{\textstyle\textstyle\sum _{i=0}^{n-(j+1)}}\left(\genfrac{}{}{0.0pt}{}{n}{j}\right)\left(\genfrac{}{}{0.0pt}{}{n-(j+1)}{i}\right){\lambda ^{i}}{\mu ^{j}}{x^{(i+j)\alpha }}{s_{k}^{n-(j+1)-i}}\Big){E_{\alpha ,\alpha }}({s_{k}})}{({s_{k}}-{s_{k+1}})({s_{k}}-{s_{k+2}})\dots ({s_{k}}-{s_{k+n}})}\hspace{3.57777pt},\hspace{3.57777pt}\hspace{3.57777pt}\hspace{3.57777pt}\]
where $k\in {\mathbb{Z}_{n+1}}$, a group under addition mod $n+1$ and, ${s_{k}}$ is a root of (12).
5 Conclusion
We have considered a random variable arising from a queueing system whose distribution is Phase-Type, and used this PH representation to derive explicit forms of the associated MML distribution. The MML class is heavy-tailed, dense in the class of all lifetime distributions on the positive real line, and mathematically tractable, making it suitable for modeling systems in which inter-arrival or service times exhibit heavy-tailed behavior. In this work, we studied the MML distributions associated with the truncated queueing systems PH/M/$1\sim n$ and M/PH/$1\sim n$, and obtained explicit density formulas when the PH distributions reduce to exponential or Erlang forms. Several examples were provided to illustrate the results. A numerical illustration was also included to demonstrate the stability of the derived expressions. The approach developed here may be extended by considering other subclasses of PH distributions, such as Coxian or generalized Erlang distributions, within the PH/M/$1\sim n$ and M/PH/$1\sim n$ frameworks. These extensions may further broaden the applicability of MML distributions in modeling complex queueing and stochastic systems.