Inter-Departure Time Correlations in PH/G/1 Queues

Abstract

In non-Markovian tandem queueing networks the output process of one site, which is the input process to the next site, is not renewal. Consequently, the correlation analysis of that output processes is essential when studying such networks. A correlation analysis in the M/G/1 queue has been studied in the literature via derivation of the joint Laplace-Stieltjes transform (LST) of the sum of two consecutive inter-departure times. That LST is obtained by considering all possible cases at departure epochs. However, those epochs are expressed via dependent variables. In this paper, we first extend the analysis to the more general PH/G/1 queue, and investigate various queues, such as E₂/G/1 and C₂/C₂/1. Then, we consider the lag-n correlation, which requires derivation of the joint LST of sum of n + 1 consecutive inter-departure times. Yet, deriving this LST by the common approach becomes impractical for n + 1 ≥ 3, as the number of all possible cases at departure epochs increases significantly. To overcome this obstacle, we derive a corresponding single-parameter LST, which expresses the sum of n + 1 consecutive inter-departure times via the (n + 1)-st departure epoch only. Consequently, the latter LST is expressed via a much fewer number of possible cases, and not less important, as a function of independent variables only, eliminating the need to derive the corresponding joint density. Considering the M/G/1 and the E₂/G/1 queues, we demonstrate that the joint LST can be reconstructed directly via the corresponding single-parameter LST when n + 1 = 2. We further conjecture that the multi-parameter joint LST can be reconstructed from the corresponding single-parameter LST in more general queues and for values of n + 1 > 2. The conjecture is validated for various PH/G/1 queues and proved for n + 1 = 3 in the M/G/1 case. The new approach facilitates the calculation of lag-n correlation of the departure process from PH/G/1 queue for n + 1 ≥ 3. Our analysis illuminates the cases when using renewal approximation of the output process provides a proper approximation when studying non-Markovian stochastic networks.

1. Introduction

Continuous-time queueing networks are rapidly growing in prevalence in modern life. Manufacturing lines, service centers, communication networks, transport systems, supply chains, etc., are just a few examples of such networks. Those networks have several connected single-site queueing systems, where jobs propagate through a network over various routes. Efficient control of such complex networks cannot be achieved without a proper probabilistic analysis. Nevertheless, except for some special cases, such as Tandem Jackson Networks [1] or Asymmetric Inclusion Process (ASIP) networks [2,3,4] that assume Markovian properties, exact expressions for the steady-state distribution of the queue lengths and performance measures of non-Markovian networks are hardly available. The Markovian property of Tandem Jackson Networks leads to a product-form solution of the joint probability-generating function of the queue lengths, as well as of the corresponding joint probability mass function. However, even under Markovian assumptions, in the ASIP model, only the joint probability-generating function of the queue workloads exhibits a product-form.

The reason that Tandem Jackson Networks have exact expressions is rooted in the fact that the departure process from each queue is also Markovian. However, in many real-life networks, inter-arrival times and service durations do not necessarily follow exponential distributions. Consequently, the inter-departure times involve correlations and the departure process becomes a complex non-renewal, requiring a further deep probabilistic analysis. Thus, the study of the departure process from single-site non-Markovian queues is primarily motivated by the quest of analyzing sequential queueing networks. Methods of queueing network analysis such as decomposition approximations usually ignore dependencies between inter-departure times and act as if the inflow from one queue to the next is a renewal process. Indeed, taking these dependencies into account has the potential for improving approximations such as in [5,6,7,8,9]. The study of the departure process from non-Markovian queues can also be applied to extend the analysis of queueing networks such as those studied in [10,11,12,13,14,15]. For example, considering the latter paper, the study of networks where products’ quality deteriorates over time can be extended to non-Markovian networks.

The steady-state probabilistic characterization of the departure instants can be expressed by the Laplace–Stieltjes transform (LST) of the inter-departure time distribution and by their correlation structure. The well-known M/G/1 queue is thoroughly studied by [16], who represents the LST of inter-departure times along with its correlation. Bitran and Dasu (1994) [17] studied the ΣPH_i/PH/1 queue and presented a formula to calculate the second moment of inter-departure times. Yeh and Chang (2000) [18] studied an M/G/1-type queue, calculated the LST of the stationary distribution of inter-departure times and provided a recursive formula for the corresponding lag-n correlation. Ferng and Chang (2001) [19] studied the lag-n correlation in a BMAP/G/1 queue. Shioda (2003) [20] derived the LST of inter-departure times, as well as their correlation, for a queue with a Markovian Arrival Process (MAP) and semi-Markovian (SM) service process. Similarly, Lim et al. (2006) [21] derived it for a queue with a Markov renewal arrival process and general service durations. Lee and Luh (2006) [22] provided the LST of inter-departure times in a PH/G/1 queue, but, due to computational complexity, restricted their numerical examples only to E_m/E_k/1 queues. Horváth et al. (2010) [23] presented the LST of inter-departure times and the stationary joint moments of two consecutive inter-departure times in an MAP/MAP/1 queue. Sagron et al. (2019) [24] provided an efficient way to express the LST of inter-departure times in a PH/G/1 queue with generalization to an ME/G/1 queue, in which the inter-arrival times follow a Matrix-Exponential (ME) distribution.

In this paper, we extend the analysis of inter-departure times from PH/G/1 queues, concentrating on the correlation aspect, as it is key for the analysis of tandem queueing networks with general service times. We focus on Phase-type (PH) distribution since it constitutes a very versatile class of distributions and can represent any ME distribution as well (see e.g., [25]). In addition, PH distributions are dense in the class of all distributions of all non-negative continuous distribution functions with any required degree of accuracy [26]. For that reason, the PH/G/1 queue covers a wide variety of systems that arise in different applications.

The correlation between two consecutive inter-departure times is derived from the joint Laplace–Stieltjes transform of the sum of two consecutive inter-departure times. The derivation of the latter joint LST for the M/G/1 queue has been already investigated in the literature by considering all possible cases at departure epochs. Those epochs are expressed via dependent variables. In this paper, we first extend this approach to derive the joint LST of two consecutive inter-departure times in the more general PH/G/1 queue. Consequently, we investigate the correlation between two consecutive inter-departure times for different levels of the squared coefficient of variation (SCV) of inter-arrival times and service duration distributions.

Then, we consider the lag-n correlation, which requires the derivation of the joint LST of the sum of n + 1 consecutive inter-departure times. Deriving this LST by the common approach becomes impractical for values of n + 1 ≥ 3, as the number of all possible cases at departure epochs increases significantly. To overcome this obstacle, we derive a corresponding single-parameter LST, which expresses the sum of n + 1 consecutive inter-departure times via the (n + 1)-st departure epoch only. This approach has two significant advantages with respect to the derivation of the latter LST: (i) it enables its derivation via much fewer possible cases, and (ii) it allows us to express it as function of independent variables only, which eliminates the need to derive the corresponding joint density, as well. We exhibit an interesting property in which the joint LST can be reconstructed directly via the corresponding single-parameter LST. This finding facilitates the derivation of the multi-parameter joint LST of the sum of n + 1 consecutive inter-departure times. For n + 1 = 2, we prove the existence of this property in M/G/1 and E₂/G/1 queues and validate it for various cases of the more general C₂/C₂/1 queue. We further conjecture that the results can be extended to the sum of n + 1 consecutive inter-departure times. To enhance this conjecture, we prove it for n + 1 = 3 in the M/G/1 queue.

Our results are directed towards a deeper analysis of non-Markovian queueing networks. In particular, calculating correlation between departure instants reveals the range in which the renewal assumption can serve as a good approximation for performance of the next site’s queue.

Also, it can possibly serve as a stepping stone to derive a better approximation in the cases where the renewal assumption fails.

The rest of the paper is as follows. For completeness, In Section 2, we give a concise description of a PH distribution. In Section 3, we derive the LST of the sum of two consecutive inter-departure times from a PH/G/1 queue, calculate the correlation between them in various queues and discuss the consequences. In Section 4, we propose a new and simpler way to express this LST for n + 1 ≥ 3 epochs. Section 5 summarizes the paper.

2. Phase-Type Distribution

For completeness, we give a short description of PH distribution as follows. A PH distribution, first introduced by [27], is the distribution of the time until absorption to state m + 1 in an (m + 1)-state continuous-time Markov chain (CTMC). In what follows, vectors and matrixes are presented by bold letters.

Let Y(t) denote the state of the CTMC at time t with the following (m + 1) $\times$ (m + 1) generator matrix Q:

$$Q=\left(\begin{array}{cc}U& U^0\\ 0& 0\end{array}\right).$$

U is an m × m nonsingular matrix with negative diagonal and non-negative off-diagonal entries and U⁰ is a non-negative column vector that satisfies U⁰ = −Ue, where e is a column vector of ones. States 1 through m are transient, such that absorption into state m + 1 from any initial state is certain. Further, let the process have an initial probability of starting in any of the m states given by the initial probability vector α = (α₁, α₂, …, α_m), where αe = 1. Namely, P(Y(0) = i) = α_i.

Let U denote the time until the process reaches the absorbing state m + 1 when starting at time t = 0 according to the m-dimensional probability vector α. That is, after starting at state i with probability α_i, the process alternates between the m transient states represented by the matrix U until finally being absorbed at state m + 1, i.e., given that Y(0) = i, U = inf{t > 0|Y(t) = m + 1}. U is said to be PH-distributed, which is denoted by PH(α,U) and called a PH representation. The cumulative distribution function (cdf) of U, F_U(u), defined for u ≥ 0 is provided by

Fu(u) = 1 − α∙exp(Uu)e.

(1)

The corresponding probability density function (pdf), f_U(u), defined for u > 0, is provided by

fu(u) = α∙exp(Uu)U⁰,

(2)

and $\tilde{U}(\theta )$, the corresponding LST defined for θ ≥ 0, is provided by

$$\tilde{U}(\theta )={\displaystyle \underset{0}{\overset{\infty }{\int }}{e}^{-\theta u}{f}_U(u)du=\mathbf{\alpha }{(\theta {\mathit{I}}_{\mathit{m}}-\mathbf{U})}^{-1}{\mathbf{U}}^0,}$$

(3)

where I_m is a unit matrix with dimension m. The nth moment of U is provided by the following (see [28]):

E(Uⁿ) = n!α(−U⁻¹)ⁿe,  n = 1, 2, …

(4)

The following theorem about PH properties, originally shown for discrete PH random variables but holding for continuous PH random variables as well, is provided in Theorem 2.6.1 in [28]:

Let Z and Y be two independent random variables. Assume that Z is PH-distributed having representation PH(α,Z) with m phases and Y is PH-distributed having representation PH(β,Y) with k phases. Then, their sum X + Y is PH-distributed having representation PH(γ,C) with m + k phases, where

$$(\gamma ,C)=\left(\left(\alpha ,\underset{k}{\underbrace{0,\hspace{0.17em}\dots \hspace{0.17em},0}}\right)\hspace{0.17em},\left(\begin{array}{cc}Z& Z^0\beta \\ 0& Y\end{array}\right)\right)$$

(5)

3. Correlation between Two Consecutive Inter-Departure Times from a PH/G/1 Queue

In this section, we analyze the correlation between two consecutive inter-departure times from a PH/G/1 queue. This is achieved by first deriving the LST of the sum of two consecutive inter-departure times and then calculating the correlation between them.

A PH/G/1 queue is a single-server infinite-buffer queue, characterized by i.i.d. PH inter-arrival times, T₁, T₂,…, all distributed as T, and by i.i.d. general service durations, X₁, X₂,…, all distributed as X. T is PH-distributed, having PH(α,T) representation with m transient phases, and LST $\tilde{T}(\theta )$; X follows a general distribution function with LST $\tilde{X}(\theta )$. It is assumed that the queue is stable, i.e., E(T) = α(−T⁻¹)e > E(X) (see [29]), namely, ρ ≡ E(X)/E(T) < 1. Let τ_n denote the departure epoch of the nth departing unit, and let D_n = τ_n − τ_n−1 denote the inter-departure time between the (n − 1)st and the nth departing unit. In addition, denote by A_n the instant when the nth unit enters service; by L_n the number of units left behind the nth departing unit; by Y_n the arrival phase of the next arriving unit at the departure instant of the nth departing unit; and by $T{r}_n^i$ the remaining inter-arrival time given Y_n = i. For ease of exposition, we count the nth departure as the 0th one.

Next, we consider a stable queue and define ${\pi }_{l,i}$ to be the probability that right after a departure instant there are l units in the system, and the next arriving unit is at its i-th phase, i = 1, …, m. The algorithm to calculate these probabilities is presented in [30].

3.1. Joint LST of the Sum of Two Consecutive Inter-Departure Times

Assume that three successive departures occur at instants τ₀, τ₁, τ₂, and consider D₁ = τ₁ − τ₀ and D₂ = τ₂ − τ₁. The sum of two inter-departure times is D₁ + D₂ = τ₂ (τ₀ = 0). The joint pdf $f_{{\tau }_2}(d_1,d_2)$ for D₁ and D₂ is denoted by (see [16])

$$f_{{\tau }_2}(d_1,d_2)\Delta d_1\Delta d_2=P(d_1\le D_1<d_1+\Delta d_1\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{d}_2\le D_2<d_2+\Delta d_2),$$

(6)

and the corresponding joint LST, ${\tilde{\tau }}_2({\theta }_1,{\theta }_2)$, is defined by

$${\tilde{\tau }}_2({\theta }_1,{\theta }_2)=E(e^{-({\theta }_1{D}_1+{\theta }_2{D}_2)})={\displaystyle \underset{d_1=0}{\overset{\infty }{\int }}{\displaystyle \underset{d_2=0}{\overset{\infty }{\int }}{e}^{-{\theta }_1{d}_1}{e}^{-{\theta }_2{d}_2}}{f}_{{\tau }_2}(d_1,}\hspace{0.17em}{d}_2)\Delta d_1\Delta d_2.$$

(7)

Given that a departure occurs at time τ₀, two successive departures can occur following one of three cases, where two of them split into two sub-cases. A departure occurring at time τ₀ leaves behind it either of the following:

(i)

No other units (L₀ = 0);

(ii)

A single unit (L₀ = 1);

(iii)

L₀ ≥ 2 units.

In cases (i) and (ii), the next departing unit, occurring at time τ₁, leads to one of two sub-cases:

(a)

It leaves no other units in the system (L₁ = 0), i.e., during its service duration no other units arrived;

(b)

It leaves L₁ ≥ 1 units behind it, i.e., during its service duration at least one unit has arrived.

Figure 1 summarizes all cases for realizations τ₁ = d₁, τ₂ = d₂.

For better clarity, we refine the number of cases into five:

(1)

Case (i), sub-case (a)—each of the two consecutive departures, occurring at times τ₀ and τ₁, leaves no other units behind it.

(2)

Case (i), sub-case (b)—a departure occurring at time τ₀ leaves no other units behind it, but a departure occurring at time τ₁ leaves L₁ ≥ 1 units behind it.

(3)

Case (ii), sub-case (a)—a departure occurring at time τ₀ leaves a single unit behind it, but a departure occurring at time τ₁ leaves no other units behind it.

(4)

Case (ii), sub-case (b)—a departure occurring at time τ₀ leaves a single unit behind it, but a departure occurring at time τ₁ leaves L₁ ≥ 1 units behind it.

(5)

Case (iii), a departure occurring at time τ₀ leaves L₀ ≥ 2 units behind it.

Thus, the sum of two inter-departure times between τ₀ and τ₂ can be represented by two instants of successive departures, as provided in Equation (8) below:

$$\begin{array}{l}{\tau }_2=D_1+D_2=\\ \\ =\left\{\begin{array}{cccc}T{r}_1^i+X_1|(T_2>X_1)+T{r}_2^j+X_2& ,& L_0=0,\hspace{0.17em}{Y}_0=i\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}(L_1=0\hspace{0.17em},\hspace{0.17em}{Y}_1=j|T_2>X_1)& i=1,\dots ,m\hspace{0.17em},\hspace{0.17em}j=1,\dots ,m\\ T{r}_1^i+X_1|(T_2\le X_1)+X_2& ,& L_0=0,\hspace{0.17em}\hspace{0.17em}{Y}_0=i,\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}(L_1\ge 1|T_2\le X_1)\hspace{0.17em}& i=1,\dots ,m\\ X_1|(T{r}_2^i>X_1)+T{r}_2^j+X_2& ,& L_0=1,\hspace{0.17em}\hspace{0.17em}{Y}_0=i\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}(L_1=0,Y_1=j|T{r}_2^j>X_1)\hspace{0.17em}& i=1,\dots ,m\hspace{0.17em},\hspace{0.17em}j=1,\dots ,m\\ X_1|(T{r}_2^i\le X_1)+X_2& ,& L_0=1,\hspace{0.17em}\hspace{0.17em}{Y}_0=i\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}(L_1\ge 1|T{r}_2^i\le X_1)& i=1,\dots ,m\\ X_1+X_2& ,& L_0\ge 2& \end{array}\right.\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\end{array}$$

(8)

In case (1), the first inter-departure time, D₁, is the sum of a remaining inter-arrival time until the next unit arrives ($T{r}_1^i$) augmented by its own service duration (X₁), given that the latter unit completes its service before the next unit arrives (T₂ > X₁). Namely, the first unit leaves no units behind it (L₁ = 0). Then, the next successive inter-departure time, D₂, is the sum of the remaining inter-arrival time until the second unit arrives ($T{r}_1^i$), augmented by its own service duration (X₂). In case (2), the first inter-departure time, D₁, is the sum of the remaining inter-arrival time until the next unit arrives ($T{r}_1^i$) augmented by its own service duration (X₁), given that the latter unit completes its service before the next unit arrives (T₂ ≤ X₁). Namely, the first unit leaves L₁ ≥ 1 units behind it. Then, the next successive inter-departure time, D₂, is the full service duration of the second departing unit (X₂). In Case (3), the first inter-departure time, D₁, is the full service duration of the first departing unit (X₁), given that it leaves no units behind it (L₁ = 0, $T{r}_2^j>X_1$). Then, the next successive inter-departure time, D₂, is the sum of the remaining inter-arrival time until the second unit arrives ($T{r}_2^j$), followed by its own service duration (X₂). In Case (4), the first inter-departure time, D₁, is the full service duration of the first departing unit (X₁), given that it leaves L₀ ≥ 1 units behind it ($T{r}_2^j\le X_1$). Then, the next successive inter-departure time, D₂, is only the service duration of the second departing unit (X₂). In Case (5), the sum of two consecutive inter-departure times, D₁ + D₂, is the sum of two service durations (X₁ + X₂). Note that the random variable τ₂ as defined in (8) is expressed in terms of 2m² + 2m + 1 random variables, not all independent. Owing to the dependencies between variables in Equation (8), we formulate the joint LST of τ₂,${\tilde{\tau }}_2({\theta }_1,{\theta }_2)$ via the joint pdf $f_{{\tau }_2}(d_1,d_2)$.

We now consider the contributions of Cases (1)–(5) to $f_{{\tau }_2}(d_1,d_2)$ as follows:

Case (1): The probability that at the departure time τ₀ = 0 the system is empty and the arrival phase of the next unit is i is provided by π_0,i. During the time until the next departure at d₁, a unit arrives between u and u + Δu with probability $f_{T{r}^i}(u)\Delta u$, so the server starts serving at time u. It completes the service between d₁ and d₁ + Δd₁, with probability f_X (d₁ − u)Δd₁. The next arriving unit starts its arrival process at time u in phase l with probability α_l. Then, the probability that during d₁ − u there are no arrivals and that the arriving unit is in its j-th arrival phase at time τ₁ = d₁ is exp(T(d₁ − u))_lj (see [28], p. 41). Similarly to d₁, the length of time until the next departure at d₂ consists of (i) the remaining time until arrival, occurring between v and v + Δv with probability $f_{T{r}^i}(v)\Delta v$, and (ii) a full service duration, starting at time v and completed between d₂ and d₂ + Δd₂ with probability f_X(d₂ − v)Δd₂. Thus, the contribution of Case (1) to $f_{{\tau }_2}(d_1,d_2)$ is

$${\displaystyle \sum _{i=1}^m{\pi }_{0,i}}{\displaystyle \sum _{l=1}^m{\displaystyle \sum _{j=1}^m{\displaystyle \underset{u=0}{\overset{d_1}{\int }}{f}_{T{r}^i}(u)f_X(d_1-u)\hspace{0.17em}{\alpha }_l\exp {(T\hspace{0.17em}(d_1-u))}_{l\hspace{0.17em}j}\Delta u}}}\hspace{0.17em}\cdot {\displaystyle \underset{v=0}{\overset{d_2}{\int }}{f}_{T{r}^j}(v)f_X(d_2-v)\Delta v}.$$

Case (2): The probability that at the departure time τ₀ = 0 the system is empty and the arrival phase of the next unit is i is provided by π_0,i. During d₁, the next unit arrives between u and u + Δu with probability $f_{T{r}^i}(u)\Delta u$, so, the server starts serving at time u. It completes the service between d₁ and d₁ + Δd₁, with probability f_X(d₁ − u)Δd₁. The probability that at least one unit arrives during the service duration d₁ − u is F_T(d₁ − u). In such a case, the server immediately starts serving the next unit at time d₁, since there are L₁ ≥ 1 units at the system. It completes the service between d₂ and d₂ + Δd₂ with probability f_X(d₂)Δd₂. Thus, the contribution of Case (2) to $f_{{\tau }_2}(d_1,d_2)$ is

$${\displaystyle \sum _{i=1}^m{\pi }_{0,i}}{\displaystyle \underset{u=0}{\overset{d_1}{\int }}{f}_{T{r}^i}(u)f_X(d_1-u)F_T(d_1-u)\Delta u}\hspace{0.17em}\cdot f_X(d_2)$$

Case (3): The probability that at the departure time τ₀ = 0 there is a single unit at the system and the next arriving unit is in its i-th arrival phase is provided by π_i,i. Since there is a single unit at the system at time τ₀, the server immediately starts serving this unit. It completes the service between d₁ and d₁ + Δd₁ with probability f_X(d₁)Δd₁. The probability that during d₁ there are no arrivals and the next unit is in its j-th arrival phase at time τ₁ = d₁ is given by exp(Td₁)_ij. During d₂, the next unit arrives between v and v + Δv with probability $f_{T{r}^i}(v)\Delta v$. The server that starts at time v completes the service between d₂ and d₂ + Δd₂, with probability f_X(d₂ − v)Δd₂. Thus, the contribution of Case (3) to $f_{{\tau }_2}(d_1,d_2)$ is

$${\displaystyle \sum _{i=1}^m{\pi }_{1,i}}{f}_X(d_1){\displaystyle \sum _{j=1}^m\exp {(T\cdot d_1)}_{i\hspace{0.17em}j}\hspace{0.17em}\cdot {\displaystyle \underset{v=0}{\overset{d_2}{\int }}{f}_{T{r}^j}(v)f_X(d_2-v)\Delta v}}$$

Case (4): The probability that at the departure time τ₀ = 0 there is a single unit at the system and the next arriving unit is in its i-th arrival phase is provided by π_1,i. The time until its arrival is T_rⁱ. Since there is a single unit at the system at time τ₀, the server immediately starts serving this unit. It completes the service between d₁ and d₁ + Δd₁ with probability f_X(d₁)Δd₁. The probability that at least one unit arrives during d₁ is $F_{T{r}^i}(d_1)$. The server starts immediately serving the next unit at time d₁, since there are L₁ ≥ 1 units at the system. It completes the service between d₂ and d₂ + Δd₂ with probability f_X(d₂)Δd₂. Thus, the contribution of case (4) to $f_{{\tau }_2}(d_1,d_2)$ is

$${\displaystyle \sum _{i=1}^m{\pi }_{1,i}}{f}_X(d_1)\hspace{0.17em}{F}_{T{r}^i}(d_1)\cdot f_X(d_2)\hspace{0.17em}$$

Case (5): The probability that at the departure time τ₀ = 0 there are at least two units in the system is provided by ${\displaystyle \sum _{l=2}^{\infty }{\displaystyle \sum _{i=1}^m{\pi }_{l,i}}}$. Then, the server starts immediately servicing two units consecutively. It completes the first service between d₁ and d₁ + Δd₁ with probability f_X(d₁)Δd₁ and the second unit between d₂ and d₂ + Δd₂ with probability f_X(d₂)Δd₂. Thus, the contribution of Case (5) to $f_{{\tau }_2}(d_1,d_2)$ is

$${\displaystyle \sum _{l=2}^{\infty }{\displaystyle \sum _{i=1}^m{\pi }_{l,i}{f}_X(d_1)\hspace{0.17em}{f}_X(d_2)\hspace{0.17em}}}$$

Collecting all the contributions of Cases (1)–(5), we obtain

$$\begin{array}{ll}{f}_{{\tau }_2}(d_1,d_2)& ={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{l=1}^m{\displaystyle \sum _{j=1}^m{\pi }_{0,i}}}}\cdot {\displaystyle \underset{u=0}{\overset{d_1}{\int }}{f}_{T{r}^i}(u)f_X(d_1-u){\alpha }_l\exp {(T(d_1-u))}_{lj}\Delta u}\hspace{0.17em}{\displaystyle \underset{v=0}{\overset{d_2}{\int }}{f}_{T{r}^j}(v)f_X(d_2-v)\Delta v}\\ & +{\displaystyle \sum _{i=1}^m{\pi }_{0,i}}{\displaystyle \underset{u=0}{\overset{d_1}{\int }}{f}_{T{r}^i}(u)f_X(d_1-u)\hspace{0.17em}{F}_T(d_1-u)\Delta u}\hspace{0.17em}\cdot f_X(d_2)\\ & +{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^m{\pi }_{1,i}}}{f}_X(d_1)\exp {(T\cdot d_1)}_{i\hspace{0.17em}j}\hspace{0.17em}\cdot {\displaystyle \underset{v=0}{\overset{d_2}{\int }}{f}_{T{r}^j}(v)f_X(d_2-v)\Delta v}\\ & +{\displaystyle \sum _{i=1}^m{\pi }_{1,i}}{f}_X(d_1)\hspace{0.17em}{F}_{T{r}^i}(d_1)\cdot f_X(d_2)\hspace{0.17em}\hspace{0.17em}+{\displaystyle \sum _{l=2}^{\infty }{\displaystyle \sum _{i=1}^m{\pi }_{l,i}}}\hspace{0.17em}{f}_X(d_1)f_X(d_2)\end{array}$$

(9)

Finally, the joint LST of the sum of two consecutive inter-departure times, ${\tilde{\tau }}_2({\theta }_1,{\theta }_2)$, is obtained by substituting Equation (9) in Equation (7).

This result generalizes Takagi’s formula for the M/G/1 queue to the PH/G/1 queue. Indeed, for exponential inter-arrival times with mean 1/λ, when substituting T = −λ and F_T (u) = $F_{T{r}^1}(u)$ = λe^−λu in Equation (9), and then substituting Equation (9) in Equation (7), Takagi’s (1991) [16] result is obtained:

$$\begin{array}{ll}{\tilde{\tau }}_2({\theta }_1,{\theta }_2)& =\hspace{0.17em}{\pi }_0\left[{\lambda }^2\frac{\tilde{X}(\lambda +{\theta }_1)}{\lambda +{\theta }_1}\frac{\tilde{X}({\theta }_2)}{\lambda +{\theta }_2}+\lambda \frac{\left[\tilde{X}({\theta }_1)-\tilde{X}({\theta }_1+\lambda )\right]}{\lambda +{\theta }_1}\tilde{X}({\theta }_2)\right]\\ & +\hspace{0.17em}\hspace{0.17em}{\pi }_1\left[\lambda \frac{\tilde{X}(\lambda +{\theta }_1)\tilde{X}({\theta }_2)}{\lambda +{\theta }_2}+\left[\tilde{X}({\theta }_1)-\tilde{X}({\theta }_1+\lambda )\right]\tilde{X}({\theta }_2)\right]\hspace{0.17em}+{\displaystyle \sum _{l=2}^{\infty }{\pi }_l}\tilde{X}({\theta }_1)\tilde{X}({\theta }_2)\hspace{0.17em}\end{array}$$

(10)

Note that in this case ${\displaystyle \sum _{i=1}^m{\pi }_{l,i}}$ shrinks to π_l.

As another example, consider the E₂/G/1 queue. Inter-arrival times are two-stage Erlang-distributed, each stage having mean 1/λ, and service durations are generally distributed. In this case, Tr¹~T~E_r(2,λ) and Tr²~Exp (λ). Substituting $T=\left(\begin{array}{cc}-\lambda & \lambda \\ 0& -\lambda \end{array}\right)$, F_T (u) = $F_{T{r}^1}(u)$ = λ²ue^−λu and $F_{T{r}^2}(u)$ = λe^−λu in (9) and then substituting (9) in (7) leads to Equation (11) below:

$$\begin{array}{l}{\tilde{\tau }}_2({\theta }_1,{\theta }_2)=\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\pi }_{0,1}\left[{\left(\frac{\lambda }{\lambda +{\theta }_1}\right)}^2\tilde{X}(\lambda +{\theta }_1)\hspace{0.17em}{\left(\frac{\lambda }{\lambda +{\theta }_2}\right)}^2-\lambda \frac{𝜕\tilde{X}(\lambda )}{𝜕\lambda }\frac{\lambda }{\lambda +{\theta }_2}+{\left(\frac{\lambda }{\lambda +{\theta }_1}\right)}^2\hspace{0.17em}\left[\tilde{X}({\theta }_1)-\tilde{X}(\lambda +{\theta }_1)+\lambda \frac{𝜕\tilde{X}(\lambda +{\theta }_1)}{𝜕(\lambda +{\theta }_1)}\right]\right]\hspace{0.17em}\tilde{X}({\theta }_2)\\ \hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}{\pi }_{0,2}\left[\frac{\lambda }{\lambda +{\theta }_1}\tilde{X}(\lambda +{\theta }_1)\hspace{0.17em}{\left(\frac{\lambda }{\lambda +{\theta }_2}\right)}^2-\frac{\lambda }{\lambda +{\theta }_1}\frac{𝜕\tilde{X}(\lambda +{\theta }_1)}{𝜕(\lambda +{\theta }_1)}\frac{\lambda }{\lambda +{\theta }_2}+\frac{\lambda }{\lambda +{\theta }_1}\hspace{0.17em}\left[\tilde{X}({\theta }_1)-\tilde{X}(\lambda +{\theta }_1)+\lambda \frac{𝜕\tilde{X}(\lambda +{\theta }_1)}{𝜕(\lambda +{\theta }_1)}\right]\right]\hspace{0.17em}\tilde{X}({\theta }_2)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\pi }_{1,1}\left[\tilde{X}(\lambda +{\theta }_1)\hspace{0.17em}{\left(\frac{\lambda }{\lambda +{\theta }_2}\right)}^2-\lambda \frac{𝜕\tilde{X}(\lambda +{\theta }_1)}{𝜕(\lambda +{\theta }_1)}\frac{\lambda }{\lambda +{\theta }_2}+\tilde{X}(\lambda )-\tilde{X}(\lambda +{\theta }_1)+\lambda \frac{𝜕\tilde{X}(\lambda +{\theta }_1)}{𝜕(\lambda +{\theta }_1)}\right]\hspace{0.17em}\tilde{X}({\theta }_2)\\ \hspace{0.17em}\hspace{0.17em}+{\pi }_{1,2}\left[\tilde{X}(\lambda )-\tilde{X}(\lambda +{\theta }_1)\right]\tilde{X}({\theta }_2)\hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{l=2}^{\infty }{\pi }_l}\tilde{X}({\theta }_1)\tilde{X}({\theta }_2)\hspace{0.17em}\end{array}$$

(11)

To further validate Equation (11), we compare it with the joint LST for various PH service distributions, estimated by simulation. In all cases studied, the analytical result is well within the confidence interval of the simulation result. Appendix A presents the simulated validation for the E₂/E₂/1 queue.

3.2. Correlations between Two Consecutive Inter-Departure Times

In this section, we analyze the correlation between two consecutive inter-departure times in the M/G/1, E₂/G/1 and C₂/G/1 queues in order to study the impact of different levels of inflow variability on the departure process.

Since E(D₁) = E(D₂) and Var(D₁) = Var(D₂), the correlation between D₁ and D₂ is provided by

$$corr(D_1,D_2)=\frac{E(D_1{D}_2)-E^2(D_1)}{\mathrm{var}(D_1)},$$

(12)

where $E(D_1{D}_2)=\frac{{𝜕}^2{\tilde{\tau }}_2({\theta }_1,{\theta }_2)}{𝜕{\theta }_1𝜕{\theta }_2}{|}_{{\theta }_1={\theta }_2=0}$ and ${\tilde{\tau }}_2({\theta }_1,{\theta }_2)$ is obtained by substituting Equation (9) in (7).

E₂ represents a two-stage Erlang distribution, having low variability, namely, SCV = 0.5 for any parameter (specifically, we chose mean 1 at each stage). C₂ represents a two-stage Coxian distribution, the variability of which depends on the distribution parameters. In order to examine the impact of high arrival variability on the departure process, we chose the following parameters: λ₁ = 4.309 for the first stage and λ₂ = 1 for the second, with p = 0.5 for the probability of moving from the first stage to the second. This implies that the arrival’s SCV in this case is SCV = 1.5. For general distribution (G), we chose the C₂ distribution with different values of SCV (less than or greater than 1). Figure 2 depicts the correlation values at each queue for the range ρ ϵ [0, 1].

Note that in the C₂/G/1 queue it is not possible to set the value of the service’s SCV within the range ρ ϵ [0, 1], since the SCV also determines the value of ρ. Therefore, in this queue, the orange line depicts a case where the SCV is greater than 1 (and not exactly 1.5) and the purple line depicts a case where the SCV is significantly greater than 1 (and not exactly 4).

Several insights can be drawn from Figure 3:

• o When the system is nearly empty (ρ→0), the departure process tends to imitate the arrival process, hence the correlation tends to zero. When the system is almost fully utilized (ρ→1), the departure process tends to imitate the service process, hence the correlation tends to zero, as well.
• o The sign of the correlation in the range ρ ϵ (0, 1) is mainly determined by the arrival process: it is negative when the arrival variability is low, and positive when it is high. This can be seen in Figure 3. In the E₂/G/1 queue with arrival’s SCV = 0.5, the correlation is negative for all examined service distributions. The negative value implies that when the inter-arrival times have relatively low variability, the departure following a short inter-departure time will be stochastically long, and vice versa. Similarly, in the C₂/G/1 queue with arrival’s SCV = 1.5, the correlation is positive for all examined service distributions. The positive value implies that when the inter-arrival times have relatively high variability, the departure following a short (long) inter-departure time will be stochastically short (long), as well.
• o In the M/G/1 queue, i.e., when the arrival’s SCV is 1, the service process has opposite impacts on the correlation sign: service distribution with SCV = 0.5 provides a positive correlation, while service distribution with SCV = 1.5 provides a negative correlation. The correlation in the case of SCV = 1 is zero since the departure from the M/M/1 queue is a renewal process.
• o Correlation analysis can help in assessing how justified the renewal assumption is as an approximation when studying the performance of queueing networks. When the correlation tends to zero, the performance calculations are more accurate, and vice versa. Interestingly, Figure 3 shows that when the service’s SCV increases, the absolute correlation value tends to zero for high utilization. Thus, the renewal assumption in these cases (high utilization with high service variability) will provide an appropriate approximation in predicting the performance of two-site tandem networks.

4. A New Approach to Obtain the Joint LST of the Sum of n + 1 Consecutive Inter-Departure Times

A higher level of correlation analysis is used by calculating the lag-n correlation, which is defined as the correlation between D₁ and D_n+1 (n = 1, 2, 3…). For this calculation, it is required to calculate the multi-parameter joint LST of the sum of n + 1 consecutive inter-departure times. It can be achieved by considering all possible scenarios expressed by dependent variables, in light of [16] and Section 3 above. However, this approach becomes very impractical as n increases. For example, when the number of phases in the arrival process is m, the corresponding joint LST, where n + 1 = 2, is expressed in terms of 2m² + 2m + 1 random variables, and where n + 1 = 3, it is expressed in terms of 2m³ + 4m² + 6m + 1 random variables.

In this section, we reveal an interesting property of the joint LST of the sum of n + 1 consecutive inter-departure times. We show that this LST can be directly reconstructed through its corresponding single-parameter LST, which is derived via much fewer possible cases and as a function of independent variables only. For example, the single-parameter LST of n + 1 = 2 is expressed in terms of 4m + 1 independent random variables, and that of n + 1 = 3 is expressed in terms of 2(4m + 1) + 2m + 1 variables. Figure 3 illustrates the number of possible cases at each approach in the cases where n + 1 = 2 and n + 1 = 3, respectively.

In this section, we reveal this property, prove its validity for some queues and validate it for various cases. First, we present an expression of the sum of two consecutive inter-departure times directly, via the departure epoch of the second departure, based only on independent variables. This leads to a direct derivation of a single-parameter LST of the aforementioned sum without calculating its joint density. Then, we show that the corresponding joint LST can be readily reconstructed through its single-parameter LST in M/G/1 and E₂/G/1 queues. We conjecture that this reconstruction property holds in more general cases. To enhance the conjecture, we validate it for various cases of the more general C₂/C₂/1 queue for n + 1 = 2, and prove it in the M/G/1 queue for n + 1 = 3. In light of this conjecture, we propose an expression of the sum of n + 1 consecutive inter-departure times directly via the (n + 1)-st departure epoch, which is based on independent variables only and leads to a direct derivation of a single-parameter LST of that sum. The reconstruction through the single-parameter LST provides the sought after multi-parameter joint LST.

4.1. Single-Parameter LST of Sum of Two Consecutive Inter-Departure Times

Although Equation (8) defines the random variable τ₂ in the PH/G/1 queue well, it is possible to express τ₂ directly via the departure epoch of the second unit and not via the sum of two consecutive inter-departure times as expressed in Equation (8). Namely, the sum of two inter-departure times can be expressed as τ₂ = A₂ + X₂, where A₂ is the length of time from τ₀ = 0 until the instant when the second unit starts service and X₂ is the latter’s service duration. This approach is much simpler than Takagi’s approach since the corresponding LST consists of fewer variables which are all independent.

The following equation defines τ₂ only as a function of the initial system states at time τ₀:

$${\tau }_2\hspace{0.17em}=\hspace{0.17em}{A}_2\hspace{0.17em}+\hspace{0.17em}{X}_2=\left\{\begin{array}{cccc}T{r}_1^i+(T_2|T_2>X_1)+X_2& ,& L_0=0,\hspace{0.17em}\hspace{0.17em}{Y}_0=i,\hspace{0.17em}\hspace{0.17em}{T}_2>X_1\hspace{0.17em}\hspace{0.17em}& i=1,\dots ,m\\ T{r}_1^i+(T_2+(X_1-T_2)|T_2\le X_1)+X_2& ,& L_0=0,\hspace{0.17em}\hspace{0.17em}{Y}_0=i,\hspace{0.17em}\hspace{0.17em}{T}_2\le X_1\hspace{0.17em}& i=1,\dots ,m\\ (T{r}_2^j|T{r}_2^j>X_1)+X_2& ,& L_0=1,\hspace{0.17em}\hspace{0.17em}{Y}_0=j,\hspace{0.17em}\hspace{0.17em}{T}_2>X_1& j=1,\dots ,m\\ (T{r}_2^j+(X_1-T{r}_2^j)|T{r}_2^j\le X_1)+X_2& ,& L_0=1,\hspace{0.17em}\hspace{0.17em}{Y}_0=j,\hspace{0.17em}\hspace{0.17em}{T}_2\le X_1\hspace{0.17em}& j=1,\dots ,m\\ X_1+X_2& ,& L_0\ge 2& \end{array}\right..$$

(13)

The explanation is as follows: When L₀ = 0, the instant A₂ is a sum of the remaining inter-arrival time until the first unit arrives, $T{r}_1^i$, plus one of two possibilities: either a full inter-arrival time of the second unit, T₂, given that the first unit completes its service before the latter unit arrives (case 1), or the sum of a full inter-arrival time of the second unit, T₂, plus (X₁ − T₂), which is the waiting time until the first unit completes its service, given that the first unit completes its service after the second unit arrives (case 2). When L₀ = 1, the instant A₂ is one of two possibilities: either a remaining inter-arrival time of the second unit, $T{r}_2^j$, given that the first unit completes its service before the next unit arrives (case 3), or the sum of a remaining inter-arrival time of the second unit, $T{r}_2^j$, plus $(X_1-T{r}_2^j)$ which is the waiting time until the first unit completes its service, given that the first unit completes its service after the second unit arrives (case 4). When L₀ ≥ 2, the instant A₂ is only X₁ (case 5).

Then, (13) becomes the following (the variables T₂ and Tr₂ in the second and forth lines in are omitted):

$${\tau }_2=\hspace{0.17em}{A}_2\hspace{0.17em}+\hspace{0.17em}{X}_2=\left\{\begin{array}{cccc}T{r}_1^i+(T_2|T_2>X_1)+X_2& ,& L_0=0,\hspace{0.17em}\hspace{0.17em}{Y}_0=i,\hspace{0.17em}\hspace{0.17em}{T}_2>X_1\hspace{0.17em}\hspace{0.17em}& i=1,\dots ,m\\ T{r}_1^i+(X_1|T_2\le X_1)+X_2& ,& L_0=0,\hspace{0.17em}\hspace{0.17em}{Y}_0=i,\hspace{0.17em}\hspace{0.17em}{T}_2\le X_1\hspace{0.17em}& i=1,\dots ,m\\ (T{r}_2^i|T{r}_2^i>X_1)+X_2& ,& L_0=1,\hspace{0.17em}\hspace{0.17em}{Y}_0=j,\hspace{0.17em}\hspace{0.17em}{T}_2>X_1& j=1,\dots ,m\\ (X_1|T{r}_2^i\le X_1)+X_2& ,& L_0=1,\hspace{0.17em}\hspace{0.17em}{Y}_0=i,\hspace{0.17em}\hspace{0.17em}{T}_2\le X_1\hspace{0.17em}& j=1,\dots ,m\\ X_1+X_2& ,& L_0\ge 2& \end{array}\right..$$

(14)

In fact, Equation (14) resembles Equation (8) but the dependencies between the variables in cases 1 and 3 in Equation (8) are replaced by independent variables in Equation (14). Furthermore, Equation (14) can be written as

$${\tau }_2=A_2+X_2=\left\{\begin{array}{cccc}T{r}_1^i+\max \hspace{0.17em}(T_2,X_1)+X_2& ,& L_0=0,\hspace{0.17em}\hspace{0.17em}{Y}_0=i& i=1,\dots ,m\\ \max \hspace{0.17em}(T{r}_2^i,X_1)+X_2& ,& L_0=1,\hspace{0.17em}\hspace{0.17em}{Y}_0=i& i=1,\dots ,m\\ X_1+X_2& ,& L_0\ge 2& \end{array}\right..$$

(15)

It should be noticed that τ₂ as defined in Equation (15) is expressed in terms of 4m + 1 independent random variables, whereas in Equation (8) it is expressed in terms of 2m² + 2m + 1 random variables, not all independent. Thus, Equation (15) facilitates a direct derivation of a single-parameter LST by a proper partitioning of possible events as follows:

$$\begin{array}{ll}{\tilde{\tau }}_2(\theta )& ={\displaystyle \sum _{i=1}^m{\pi }_{0,i}}\left[E(e^{-\theta T{r}_1^i})\left(E(e^{-\theta \hspace{0.17em}{T}_2};T_2>X_1)+E(e^{-\theta X_1};T_2\le X_1)\right)E(e^{-\theta X_2})\right]\\ & +\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{j=1}^m{\pi }_{1,j}}\left[\left(E(e^{-\theta \hspace{0.17em}{T}_2};T_2>X_1)+E(e^{-\theta X_1};T_2\le X_1)\right)E(e^{-\theta X_2})\right]+{\displaystyle \sum _{l=2}^{\infty }{\displaystyle \sum _{i=1}^m{\pi }_{l,i}}}E(e^{-\theta X_1})E(e^{-\theta X_2})\end{array},$$

(16)

In other words, (16) is written as

$$\begin{array}{ll}{\tilde{\tau }}_2(\theta )& ={\displaystyle \sum _{i=1}^m{\pi }_{0,i}}\left[E(e^{-\theta T{r}_1^i})\left({\displaystyle \underset{0}{\overset{\infty }{\int }}\left(f_{X_1}(x){\displaystyle \underset{x}{\overset{\infty }{\int }}{e}^{-\theta \cdot t}{f}_{T_2}(t)dt}\right)dx}+{\displaystyle \underset{0}{\overset{\infty }{\int }}\left(f_{T_2}(t){\displaystyle \underset{t}{\overset{\infty }{\int }}{e}^{-\theta \cdot x}{f}_{X_1}(x)dx}\right)dt}\right)E(e^{-\theta X_2})\right]\\ & +\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{j=1}^m{\pi }_{1,j}}\left[\left({\displaystyle \underset{0}{\overset{\infty }{\int }}\left(f_{X_1}(x){\displaystyle \underset{x}{\overset{\infty }{\int }}{e}^{-\theta \cdot t}{f}_{T{r}_2^j}(t)dt}\right)dx}+{\displaystyle \underset{0}{\overset{\infty }{\int }}\left(f_{T{r}_2^j}(t){\displaystyle \underset{t}{\overset{\infty }{\int }}{e}^{-\theta \cdot x}{f}_{X_1}(x)dx}\right)dt}\right)E(e^{-\theta X_2})\right]\\ & +{\displaystyle \sum _{l=2}^{\infty }{\displaystyle \sum _{i=1}^m{\pi }_{l,i}}}E(e^{-\theta X_1})E(e^{-\theta X_2})\end{array}.$$

(17)

In the case of M/G/1 queue, Equation (17) becomes

$$\begin{array}{ll}{\tilde{\tau }}_2(\theta )& =\hspace{0.17em}{\pi }_0\left[{\tilde{T}}_1(\theta )\left({\tilde{T}}_2(\theta ){\tilde{X}}_1(\lambda +\theta )+{\tilde{X}}_1(\theta )-{\tilde{X}}_1(\lambda +\theta )\right){\tilde{X}}_2(\theta )\right]\\ & +\hspace{0.17em}\hspace{0.17em}{\pi }_1\hspace{0.17em}\left[\left({\tilde{T}}_2(\theta ){\tilde{X}}_1(\lambda +\theta )+{\tilde{X}}_1(\theta )-{\tilde{X}}_1(\lambda +\theta )\right){\tilde{X}}_2(\theta )\hspace{0.17em}\right]+{\displaystyle \sum _{l=2}^{\infty }{\pi }_l}{\tilde{X}}_1(\theta ){\tilde{X}}_2(\theta )\end{array}.$$

(18)

When θ₁ = θ₂ = θ, the single-parameter LST of the sum of two consecutive inter-departure times, Equation (18), coincides with the joint LST of the sum of two consecutive inter-departure times, Equation (10)—Takagi’s (1991) [16] result.

In the case of E₂/G/1 queue, Equation (17) becomes

$$\begin{array}{l}{\tilde{\tau }}_2(\theta )=\hspace{0.17em}\\ =\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\pi }_{0,1}\left[{\tilde{T}}_1(\theta )\tilde{X}(\lambda +\theta ){\tilde{T}}_2(\theta )-\lambda \frac{𝜕\tilde{X}(\lambda )}{𝜕\lambda }\tilde{T}{r}_1(\theta )+\tilde{T}{r}_1(\theta )\left[\tilde{X}(\theta )-\tilde{X}(\lambda +\theta )+\lambda \frac{𝜕\tilde{X}(\lambda +\theta )}{𝜕(\lambda +\theta )}\right]\right]{\tilde{X}}_2(\theta )\\ \hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}{\pi }_{0,2}\left[\tilde{T}{r}_1(\theta )\tilde{X}(\lambda +\theta ){\tilde{T}}_2(\theta )-\tilde{T}{r}_1(\theta )\frac{𝜕\tilde{X}(\lambda +\theta )}{𝜕(\lambda +\theta )}\tilde{T}{r}_2(\theta )+\tilde{T}{r}_1(\theta )\left[\tilde{X}(\theta )-\tilde{X}(\lambda +\theta )+\lambda \frac{𝜕\tilde{X}(\lambda +\theta )}{𝜕(\lambda +\theta )}\right]\right]{\tilde{X}}_2(\theta )\\ \hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\pi }_{1,1}\left[\tilde{X}(\lambda +\theta ){\tilde{T}}_2(\theta )-\lambda \frac{𝜕\tilde{X}(\lambda +\theta )}{𝜕(\lambda +\theta )}\tilde{T}{r}_2(\theta )+\tilde{X}(\lambda )-\tilde{X}(\lambda +\theta )+\lambda \frac{𝜕\tilde{X}(\lambda +\theta )}{𝜕(\lambda +\theta )}\right]{\tilde{X}}_2(\theta )\\ \hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}{\pi }_{1,2}\left[\tilde{X}(\lambda )-\tilde{X}(\lambda +\theta )\right]{\tilde{X}}_2(\theta )\hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{l=2}^{\infty }{\pi }_l}{\tilde{X}}_1(\theta ){\tilde{X}}_2(\theta )\end{array}.$$

(19)

The importance of (18) and (19) is that the joint LSTs—(10) and (11), respectively—can be reconstructed by substituting θ_i in (18) and in (19) instead of θ, as follows: ${\tilde{T}}_i(\theta )$ is replaced by ${\tilde{T}}_i({\theta }_i)$, $\tilde{T}{r}_i(\theta )$ is replaced by $\tilde{T}{r}_i({\theta }_i)$, ${\tilde{X}}_i(\theta )$ is replaced by ${\tilde{X}}_i({\theta }_i)$ and $X_i(\lambda +\theta )$ is replaced by ${\tilde{X}}_i(\lambda +{\theta }_i)$. Thus, the joint LST of the sum of two consecutive inter-departure times in M/G/1 queue can be constructed from (18) as follows:

$$\begin{array}{l}{\tilde{\tau }}_2({\theta }_1,{\theta }_2)=\\ =\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\pi }_0\left[{\tilde{T}}_1({\theta }_1)\left({\tilde{T}}_2({\theta }_2){\tilde{X}}_1(\lambda +{\theta }_1)+{\tilde{X}}_1({\theta }_1)-{\tilde{X}}_1(\lambda +{\theta }_1)\right){\tilde{X}}_2({\theta }_2)\right]\\ \hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}{\pi }_1\hspace{0.17em}\left[\left({\tilde{T}}_2({\theta }_2){\tilde{X}}_1(\lambda +{\theta }_1)+{\tilde{X}}_1({\theta }_1)-{\tilde{X}}_1(\lambda +{\theta }_1)\right){\tilde{X}}_2({\theta }_2)\hspace{0.17em}\right]+{\displaystyle \sum _{l=2}^{\infty }{\pi }_l}\hspace{0.17em}{\tilde{X}}_1({\theta }_1){\tilde{X}}_2({\theta }_2)\hspace{0.17em}\end{array},$$

(20)

which coincides with Equation (10) (note that ${\tilde{T}}_i({\theta }_i)=\frac{\lambda }{\lambda +{\theta }_i}$; ${\tilde{X}}_i(\lambda +{\theta }_1)=\tilde{X}(\lambda +{\theta }_1)$; and ${\tilde{X}}_i({\theta }_i)=\tilde{X}({\theta }_i)$ for i = 1, 2). Similarly, the joint LST of the sum of two consecutive inter-departure times in the E₂/G/1 queue can be constructed from (18) as follows:

$$\begin{array}{l}{\tilde{\tau }}_2({\theta }_1,{\theta }_2)=\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\pi }_{0,1}\left[{\tilde{T}}_1({\theta }_1)\tilde{X}(\lambda +{\theta }_1){\tilde{T}}_2({\theta }_2)-\lambda \frac{𝜕\tilde{X}(\lambda )}{𝜕\lambda }\tilde{T}{r}_1({\theta }_1)+\tilde{T}{r}_1({\theta }_1)\left[\tilde{X}({\theta }_1)-\tilde{X}(\lambda +{\theta }_1)+\lambda \frac{𝜕\tilde{X}(\lambda +{\theta }_1)}{𝜕(\lambda +{\theta }_1)}\right]\right]{\tilde{X}}_2({\theta }_2)\\ \hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}{\pi }_{0,2}\left[\tilde{T}{r}_1({\theta }_1)\tilde{X}(\lambda +{\theta }_1){\tilde{T}}_2({\theta }_2)-\tilde{T}{r}_1({\theta }_1)\frac{𝜕\tilde{X}(\lambda +{\theta }_1)}{𝜕(\lambda +{\theta }_1)}\tilde{T}{r}_2({\theta }_2)+\tilde{T}{r}_1({\theta }_1)\left[\tilde{X}({\theta }_1)-\tilde{X}(\lambda +{\theta }_1)+\lambda \frac{𝜕\tilde{X}(\lambda +{\theta }_1)}{𝜕(\lambda +{\theta }_1)}\right]\right]{\tilde{X}}_2({\theta }_2)\\ \hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}{\pi }_{1,1}\left[\tilde{X}(\lambda +{\theta }_1){\tilde{T}}_2({\theta }_2)-\lambda \frac{𝜕\tilde{X}(\lambda +{\theta }_1)}{𝜕(\lambda +{\theta }_1)}\tilde{T}{r}_2({\theta }_2)+\tilde{X}(\lambda )-\tilde{X}(\lambda +{\theta }_1)+\lambda \frac{𝜕\tilde{X}(\lambda +{\theta }_1)}{𝜕(\lambda +{\theta }_1)}\right]{\tilde{X}}_2({\theta }_2)\\ \hspace{0.17em}+{\pi }_{1,2}\left[\tilde{X}(\lambda )-\tilde{X}(\lambda +{\theta }_1)\right]{\tilde{X}}_2({\theta }_2)\hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{l=2}^{\infty }{\pi }_l}{\tilde{X}}_1({\theta }_1){\tilde{X}}_2({\theta }_2)\end{array},$$

(21)

which coincides with Equation (11) (here, ${\tilde{T}}_i({\theta }_i)={\left(\frac{\lambda }{\lambda +{\theta }_i}\right)}^2$; $\tilde{T}{r}_i({\theta }_i)=\frac{\lambda }{\lambda +{\theta }_i}$; ${\tilde{X}}_i(\lambda +{\theta }_1)=\tilde{X}(\lambda +{\theta }_1)$; and ${\tilde{X}}_i({\theta }_i)=\tilde{X}({\theta }_i)$ for i = 1, 2).

We conjecture that the same reconstruction can be used in other PH/G/1 queues, as well. For that purpose, we validate this reconstruction in the case of the C₂/C₂/1 queue for various parameters, as the Coxian distribution can represent any PH distribution [31] with much fewer parameters [32]. We compare the reconstructed joint LST obtained by Equation (19) vs. the joint LST estimated by simulation. In all cases examined, the reconstructed joint LST is well within the confidence interval of the joint LST estimated by simulation. Figure 4 presents the case of a C₂/C₂/1 queue where the inter-arrival times and service durations are two-stage Coxian-distributed. The parameters of inter-arrival times are as follows: λ₁ = 1 for the first stage and λ₂ = 2 for the second, with p = 1 for the probability of moving from the first stage to the second. The parameters of service durations are as follows: μ₁ = 1 for the first stage and μ₂ = 2 for the second, with p = 1 for the probability of moving from the first stage to the second. Each simulation run is 10⁶ departing units after a warm-up period of 5 × 10⁴ units. For a given θ₁ and θ₂, the joint LST is estimated by

$${\tilde{\widehat{\tau }}}_2({\theta }_1,{\theta }_2)=\frac{{\displaystyle \sum _{i=1}^{{10}^6}{e}^{-({\theta }_1{\tau }_i+{\theta }_2{\tau }_{i+!})}}}{{10}^6},$$

(22)

where τ_i is an i-th inter-departure time obtained by simulation. To show the results in a two-dimensional graphical exposition, we exhibit two examples: (i) θ₁= θ₂ = θ and (ii) θ₁ = θ, θ₂ = 0. Note that the latter is exactly the marginal LST of τ₁. The corresponding graphs are depicted in Figure 4.

4.2. Sum of n + 1 Consecutive Inter-Departure Times in a PH/G/1 Queue

The LST of the sum of n + 1 consecutive inter-departure times is constructed via the joint pdf, $f_{{\tau }_{n+1}}(d_1,\dots ,d_{n+1})$ for D₁, …, D_{n + 1}, which is defined by

$$f_{{\tau }_{n+1}}(d_1,\dots ,d_{n+1})\Delta d_1\dots \Delta d_{n+1}=P(d_1<D_1<d_1+\Delta d_1\hspace{0.17em},\dots ,\hspace{0.17em}{d}_{n+1}<D_{n+1}<d_{n+1}+\Delta d_{n+1}).$$

(23)

Then, the corresponding joint LST is

$${\tilde{\tau }}_{n+1}({\theta }_1,\dots ,{\theta }_{n+1})=E(e^{-{\displaystyle \sum _{i=1}^{n+1}{\theta }_i{D}_i}})={\displaystyle \underset{d_1=0}{\overset{\infty }{\int }}\dots {\displaystyle \underset{d_{n+1}=0}{\overset{\infty }{\int }}{e}^{-{\displaystyle \sum _{i=1}^{n+1}{\theta }_i{d}_i}}}{f}_{{\tau }_n+1}\hspace{0.17em}{d}_{n+1})\Delta d_1\dots \Delta d_{n+1}}\hspace{0.17em},$$

(24)

The joint pdf, Equation (22), can be obtained by extending Equation (9) to the case of n + 1 ≥ 3, but the amount of possible cases which involve dependent variables swells as n increases as demonstrated in Figure 3 and in Appendix B. Therefore, this procedure becomes impractical. However, the sum of n + 1 consecutive inter departure times, τ_n+1, can be expressed directly via the departure epoch of the (n + 1)-st departing unit (n + 1 > 2), as follows:

$${\tau }_{n+1}=\hspace{0.17em}{A}_{n+1}\hspace{0.17em}+\hspace{0.17em}{X}_{n+1}=\left\{\begin{array}{cccc}T{r}_1^i+\max \hspace{0.17em}\hspace{0.17em}\left({\displaystyle \sum _{l=2}^{n+1}{T}_l},A_{i,0,n}\right)+X_{n+1}& ,& L_0=0\hspace{0.17em},\hspace{0.17em}{Y}_0=i& i=1,\dots ,m\\ ⋮& ,& ⋮& ⋮\\ \max \hspace{0.17em}\hspace{0.17em}\left(T{r}_{j+1}^i+{\displaystyle \sum _{l=j+2}^{n+1}{T}_l},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{A}_{i,j,n}\right)+X_{n+1}& ,& L_0=j\hspace{0.17em},\hspace{0.17em}{Y}_0=i& i=1,\dots ,m,\hspace{0.17em}\hspace{0.17em}j=1,\dots ,n\\ ⋮& ,& ⋮& ⋮\\ {\displaystyle \sum _{l=1}^{n+1}{X}_l}& ,& L_0\ge n+1& \end{array}\right.,$$

(25)

where A_i,j,n, representing the length of time from τ₀ = 0 until the instant when the nth unit starts service, is defined recursively (A_i,j,1 = 0, X₀ = 0):

$$A_{i,j,n}=\left\{\begin{array}{ccc}\max \hspace{0.17em}\left({\displaystyle \sum _{l=2}^n{T}_l},A_{i,j,n-1}\right)+X_n& ,& j=0\\ \max \hspace{0.17em}\hspace{0.17em}\left(T{r}_{j+1}^i+{\displaystyle \sum _{l=j+2}^n{T}_l},A_{i,j,n-1}\right)+X_n& ,& 0<j<n\\ {\displaystyle \sum _{l=1}^{n-1}{X}_l}& ,& j=n\end{array}\right..$$

(26)

Thus, ${\tilde{\tau }}_n(\theta )$ can be obtained directly by a proper partitioning of all possible events as a function of independent variables, as demonstrated in Equation (17). We note that ${\tilde{\tau }}_n(\theta )$ holds for any service time X, while the computational effort increases with n.

To derive the corresponding joint LST, ${\tilde{\tau }}_n({\theta }_1,\dots ,{\theta }_n)$, we conjecture the following, in light of the results presented in Section 4.1.

Conjecture 1.

The joint LST of the sum of n + 1 consecutive inter-departure times in the PH/G/1 queue,${\tilde{\tau }}_n({\theta }_1,\dots ,{\theta }_n)$, can be reconstructed via its single-parameter LST, ${\tilde{\tau }}_n(\theta )$, by substituting in ${\tilde{\tau }}_n(\theta )$: (i) ${\tilde{T}}_i({\theta }_i)$ instead of ${\tilde{T}}_i(\theta )$, and (ii) ${\tilde{X}}_i({\theta }_i)$ instead of ${\tilde{X}}_i(\theta )$, i = 1, 2,…, n.

For example, consider the case n + 1 = 3. First, we present τ₃ directly via the departure epoch of the third departing unit, so the single-parameter LST is expressed in terms of 2(4m + 1) + 2m + 1 independent random variables. Then, we obtain its joint LST with no need for the derivation of the three-dimensional joint pdf. The proof of the conjecture for the joint LST of n + 1 = 3 consecutive inter-departure times in the M/G/1 queue appears in Appendix C.

From Equation (25) we express τ_{3 =} A₃ + X₃ as follows:

$${\tau }_3\hspace{0.17em}=A_3\hspace{0.17em}+X_3=\left\{\begin{array}{cccc}T{r}_1^i+\max (T_2+T_3,\max (T_2,X_1)+X_2)+X_3& ,& L_0=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{Y}_0=i& i=1,\dots ,m\\ \max (T{r}_2^i+T_3,\max (T{r}_2^i,X_1)+X_2)+X_3& ,& L_0=1,\hspace{0.17em}{Y}_0=i& i=1,\dots ,m\\ \max (T{r}_3^i,X_1+X_2)+X_3& ,& L_0=2,\hspace{0.17em}{Y}_0=i& i=1,\dots ,m\\ X_1+X_2+X_3& ,& L_0\ge 3& \end{array}\right..$$

(27)

This representation depends only on the initial states, which, by proper partitioning of possible events, eliminates the need to derive the joint pdf and facilitates a direct derivation of a single-parameter LST as follows:

$$\begin{array}{l}{\tilde{\tau }}_3(\theta )=\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{i=1}^m{\pi }_{0,i}}\left[E(e^{-\theta \hspace{0.17em}\hspace{0.17em}T{r}_1^i})\left(\begin{array}{l}E(e^{-\theta \hspace{0.17em}\hspace{0.17em}{T}_2};T_2>X_1)\left(E(e^{-\theta \hspace{0.17em}{X}_2};T_3\le X_2)+E(e^{-\theta \hspace{0.17em}{T}_3};T_3>X_2)\right)+\\ E(e^{-\theta (\hspace{0.17em}{X}_1+X_2)};T_3\le X_2,T_2\le X_1)+E(e^{-\theta (\hspace{0.17em}{T}_2+T_3)};T_3\le X_2,T_2\le X_1,T_2+T_3>X_1+X_2)+\\ E(e^{-\theta (\hspace{0.17em}{X}_1+X_2)};T_3\le X_2,T_2\le X_1,T_2+T_3\le X_1+X_2)\end{array}\right)E(e^{-\theta X_3})\right]\\ +\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{i=1}^m{\pi }_{1,i}}\left[\left(\begin{array}{l}E(e^{-\theta \hspace{0.17em}\hspace{0.17em}T{r}_2^i};T{r}_2^i>X_1)\left(E(e^{-\theta (\hspace{0.17em}{X}_2};T_3\le X_2)+E(e^{-\theta (T_3};T_3>X_2)\right)+\\ E(e^{-\theta (\hspace{0.17em}{X}_1+X_2)};T_3\le X_2,T{r}_2^i\le X_1)+E(e^{-\theta (\hspace{0.17em}{T}_2+T_3)};T_3\le X_2,T_2>X_1,T{r}_2^i+T_3>X_1+X_2)\\ E(e^{-\theta (\hspace{0.17em}{X}_1+X_2)};T_3\le X_2,T_2>X_1,T{r}_2^i+T_3\le X_1+X_2)\end{array}\right)E(e^{-\theta X_3})\right]\hspace{0.17em}\hspace{0.17em}\\ +\hspace{0.17em}\hspace{0.17em}{\displaystyle \sum _{i=1}^m{\pi }_{2,i}\left[\left(E(e^{-\theta \hspace{0.17em}\hspace{0.17em}T{r}_3^i};T{r}_3^i>X_1+X_2)+E(e^{-\theta \hspace{0.17em}(X_1+X_2)};T{r}_3^i\le X_1+X_2)\right)E(e^{-\theta X_3})\right]}\hspace{0.17em}+{\displaystyle \sum _{l=2}^{\infty }{\displaystyle \sum _{i=1}^m{\pi }_{l,i}}}E(e^{-\theta X_1})E(e^{-\theta X_2})E(e^{-\theta X_3})\end{array}.$$

(28)

By applying the conjecture, one can reconstruct the joint LST of three consecutive inter-departure times via Equation (28).

For example, in case of the M/G/1 queue, Equation (28) is reduced to

$$\begin{array}{l}{\tilde{\tau }}_3(\theta )=\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\pi }_0\left[{\tilde{T}}_1({\theta }_1)\left(\begin{array}{l}{\tilde{T}}_2(\theta ){\tilde{X}}_1(\lambda +\theta )\left({\tilde{X}}_2(\theta )-{\tilde{X}}_2(\lambda +\theta )+{\tilde{T}}_3(\theta ){\tilde{X}}_2(\lambda +\theta )\right)+\\ \left({\tilde{X}}_1(\theta )-{\tilde{X}}_1(\lambda +\theta )\right)\left({\tilde{X}}_2(\theta )-{\tilde{X}}_2(\lambda +\theta )\right)+\\ {\tilde{T}}_3(\theta )\lambda {\tilde{X}}_2(\lambda +\theta )\frac{𝜕{\tilde{X}}_1(\lambda +\theta )}{𝜕(\lambda +\theta )}+\\ {\tilde{X}}_1(\theta ){\tilde{X}}_2(\lambda +\theta )-{\tilde{X}}_1(\lambda +\theta ){\tilde{X}}_2(\lambda +\theta )-\lambda {\tilde{X}}_2(\lambda +\theta )\frac{𝜕{\tilde{X}}_1(\lambda +\theta )}{𝜕(\lambda +\theta )}\end{array}\right){\tilde{X}}_3(\theta )\right]\\ \hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}{\pi }_1\left[\left(\begin{array}{l}{\tilde{T}}_2(\theta ){\tilde{X}}_1(\lambda +\theta )\left({\tilde{X}}_2(\theta )-{\tilde{X}}_2(\lambda +\theta )+{\tilde{T}}_3(\theta ){\tilde{X}}_2(\lambda +\theta )\right)+\\ \left({\tilde{X}}_1(\theta )-{\tilde{X}}_1(\lambda +\theta )\right)\left(\tilde{X}(\theta )-{\tilde{X}}_2(\lambda +\theta )\right)+\\ {\tilde{T}}_3(\theta )\lambda {\tilde{X}}_2(\lambda +\theta )\frac{𝜕{\tilde{X}}_1(\lambda +\theta )}{𝜕(\lambda +\theta )}+\\ {\tilde{X}}_1(\theta ){\tilde{X}}_2(\lambda +\theta )-{\tilde{X}}_1(\lambda +\theta ){\tilde{X}}_2(\lambda +\theta )-\lambda {\tilde{X}}_2(\lambda +\theta )\frac{𝜕{\tilde{X}}_1(\lambda +\theta )}{𝜕(\lambda +\theta )}\end{array}\right){\tilde{X}}_3(\theta )\right]\\ \hspace{0.17em}\hspace{0.17em}+{\pi }_2\left[\left({\tilde{T}}_3(\theta ){\tilde{X}}_1(\lambda +\theta ){\tilde{X}}_2(\lambda +\theta )+\left({\tilde{X}}_1(\theta )-{\tilde{X}}_1(\lambda +\theta )\right)\left({\tilde{X}}_2(\theta )-{\tilde{X}}_2(\lambda +\theta )\right)\right){\tilde{X}}_3(\theta )\right]\hspace{0.17em}\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}+{\displaystyle \sum _{l=3}^{\infty }{\pi }_l}{\tilde{X}}_1(\theta ){\tilde{X}}_2(\theta ){\tilde{X}}_3(\theta )\end{array}.$$

(29)

According to the Conjecture, the joint LST, ${\tilde{\tau }}_3({\theta }_1,{\theta }_2,{\theta }_3)$, can be reconstructed via Equation (29), by substituting the following: (i) ${\tilde{T}}_i({\theta }_i)$ instead of ${\tilde{T}}_i(\theta )$; (ii) ${\tilde{X}}_i({\theta }_i)$ instead of ${\tilde{X}}_i(\theta )$; and (iii) ${\tilde{X}}_1(\lambda +{\theta }_1)$ instead of ${\tilde{X}}_1(\lambda +\theta )$. Thus, the joint LST ${\tilde{\tau }}_3({\theta }_1,{\theta }_2,{\theta }_3)$ is derived from Equation (21) as follows:

$$\begin{array}{l}{\tilde{\tau }}_3({\theta }_1,{\theta }_2,{\theta }_3)=\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\pi }_0\left[{\tilde{T}}_1({\theta }_1)\left(\begin{array}{l}{\tilde{T}}_2({\theta }_1){\tilde{X}}_1(\lambda +{\theta }_1)\left({\tilde{X}}_2({\theta }_2)-{\tilde{X}}_2(\lambda +{\theta }_2)+{\tilde{T}}_3({\theta }_3){\tilde{X}}_2(\lambda +{\theta }_2)\right)+\\ \left({\tilde{X}}_1({\theta }_1)-{\tilde{X}}_1(\lambda +{\theta }_1)\right)\left({\tilde{X}}_2({\theta }_2)-{\tilde{X}}_2(\lambda +{\theta }_2)\right)+\\ {\tilde{T}}_3({\theta }_3)\lambda {\tilde{X}}_2(\lambda +{\theta }_2)\frac{𝜕{\tilde{X}}_1(\lambda +{\theta }_1)}{𝜕(\lambda +{\theta }_1)}+\\ {\tilde{X}}_1({\theta }_1){\tilde{X}}_2(\lambda +{\theta }_2)-{\tilde{X}}_1(\lambda +{\theta }_1){\tilde{X}}_2(\lambda +{\theta }_2)-\lambda {\tilde{X}}_2(\lambda +{\theta }_2)\frac{𝜕{\tilde{X}}_1(\lambda +{\theta }_1)}{𝜕(\lambda +{\theta }_1)}\end{array}\right){\tilde{X}}_3({\theta }_3)\right]\\ \hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}{\pi }_1\left[\left(\begin{array}{l}{\tilde{T}}_2({\theta }_2){\tilde{X}}_1(\lambda +{\theta }_1)\left({\tilde{X}}_2({\theta }_2)-{\tilde{X}}_2(\lambda +{\theta }_2)+{\tilde{T}}_3({\theta }_3){\tilde{X}}_2(\lambda +{\theta }_2)\right)+\\ \left({\tilde{X}}_1({\theta }_1)-{\tilde{X}}_1(\lambda +{\theta }_1)\right)\left({\tilde{X}}_2(\theta )-{\tilde{X}}_2(\lambda +{\theta }_2)\right)+\\ {\tilde{T}}_3({\theta }_3)\lambda {\tilde{X}}_2(\lambda +{\theta }_2)\frac{𝜕{\tilde{X}}_1(\lambda +{\theta }_1)}{𝜕(\lambda +{\theta }_1)}+\\ {\tilde{X}}_1({\theta }_1){\tilde{X}}_2(\lambda +{\theta }_2)-{\tilde{X}}_1(\lambda +{\theta }_1){\tilde{X}}_2(\lambda +{\theta }_2)-\lambda {\tilde{X}}_2(\lambda +{\theta }_2)\frac{𝜕{\tilde{X}}_1(\lambda +{\theta }_1)}{𝜕(\lambda +{\theta }_1)}\end{array}\right){\tilde{X}}_3({\theta }_3)\right]\\ \hspace{0.17em}\hspace{0.17em}+{\pi }_2\left[\left({\tilde{T}}_3({\theta }_3){\tilde{X}}_1(\lambda +{\theta }_1){\tilde{X}}_2(\lambda +{\theta }_2)+\left({\tilde{X}}_1({\theta }_1)-{\tilde{X}}_1(\lambda +{\theta }_1)\right)\left({\tilde{X}}_2({\theta }_2)-{\tilde{X}}_2(\lambda +{\theta }_2)\right)\right){\tilde{X}}_3({\theta }_3)\right]\hspace{0.17em}\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}+{\displaystyle \sum _{l=3}^{\infty }{\pi }_l}{\tilde{X}}_1({\theta }_1){\tilde{X}}_2({\theta }_2){\tilde{X}}_3({\theta }_3)\end{array}.$$

(30)

The proof of the conjecture in this case appears in Appendix C.

5. Summary

This paper concentrates on the analysis of inter-departure time correlations from the PH/G/1 queue. This analysis is aimed at enabling a more accurate investigation of the performance of non-Markovian tandem queues, where the output process from one site is the input process to the next. We first derive the joint LST of the sum of two inter-departure times of the PH/G/1 queue by considering all possible cases at departure epochs in light of [16] approach for the M/G/1 queue. Consequently, the correlation between two consecutive inter-departure times for various queues is calculated and investigated. This investigation indicates cases when using the renewal assumption of the output process provides a proper approximation when studying the performance of tandem queueing networks. Since the derivation of the joint LST by the abovementioned common approach becomes impractical for n + 1 ≥ 3, a new approach is proposed by which the multi-parameter joint LST of the sum of n + 1 consecutive inter-departure times is reconstructed via the corresponding single-parameter LST. This approach is based on much fewer possible cases, all expressed in terms of independent variables. Consequently, the lag-n correlation can be calculated and used for the investigation of performance assessment in non-Markovian tandem queueing networks.