A dual tandem queueing system with GI service time at the first queue

In this paper we consider the analysis of a tandem queueing model 
$M/G/1 -> ./M/1$. In contrast to the vast majority of the previous 
literature on tandem queuing models we consider the case with GI 
service time at the first queue and with infinite buffers. The 
system can be described by an M/G/1-type Markov process at the 
departure epochs of the first queue. The main result of the paper is 
the steady-state vector generating function at the embedded epochs, 
which characterizes the joint distribution of the number of 
customers at both queues. The steady-state Laplace-Stieljes 
transform and the mean of the sojourn time of the customers in the 
system are also obtained. 
   
We provide numerical examples and discuss the dependency of the 
steady-state mean of the sojourn time of the customers on several 
basic system parameters. Utilizing the structural characteristics of 
the model we discuss the interpretation of the results. This gives 
an insight into the behavior of this tandem queuing model and can be 
a base for developing approximations for it.


1.
Introduction. Tandem queueing model consists of sequentially interconnected queues and it is a special case of open queueing networks. This model became great attention in the last decades, since it has several application areas. For example it is commonly used in telecommunication networks (see e.g. in [1]).
Tandem queueing models with infinite and finite buffers has huge literature. The vast majority of the literature deals with the case of Markovian service times, like e.g. the work [6] in which a tandem queuing model with Discrete-Time Markovian Arrival Process (D-MAP) and Phase Type (PH) service times is investigated. Another example is [8], which deals with a multi-server model with finite number of queues, MAP arrival and exponential service times. This model has no buffer at all. Numerous approximation methods have been developed for tandem queuing model with Markovian service times. One method is based on the decomposition of the system into single-server finite buffer stations ( [5], [16]). Another approach is the product-form approximation of model with infinite buffers ( [3]). Tandem queueing model with deterministic service times is analyzed by Boxma and Resing in [2]. A very general model with GI service time at the first server and with finite intermediate buffer is investigated and proposed as model of a call center by Dudin et al. in [4].
All of these works consider either specific (Markovian or deterministic) service time at the first queue or at least one finite buffer in the model.
In contrast to the above references, in this paper we consider the case with GI service time at the first queue and with infinite buffers. Analysis of infinite buffer system is somewhat simpler from theoretical point of view comparing to the finite buffer counterpart of the model. This rise the chance to get at least partly closedform results, whose interpretation can contribute to the deeper understanding of the behavior of the model. The considered system can be described by an M/G/1type Markov chain embedded at the customer departure epochs of the first queue. We apply standard probability-generating function (PGF) techniques and matrix analytic methods for the analysis of the system.
The contribution of this paper is the analysis of the tandem queueing model M/G/1− > ./M/1. We derive the governing equation of the system in terms of vector generating functions (vector GF) and matrix generating functions. This relation characterizes the evolution of the system at the embedded epochs. The main result of the paper is the steady-state vector GF, which characterizes the joint distribution of the the number of customers at both stations. Based on it we obtain also the steady-state probability vector of the number of customers at the second station and the vector factorial moments of the steady-state number of customers at the first station, both at the embedded epochs. The steady-state Laplace-Stieljes transform (LST) and the steady-state mean of the sojourn time of the customers in the system are also given.
We discuss the computational procedure and we provide numerical examples to demonstrate it and to have a short look into the dependency of the steady-state mean of the sojourn time of the customers on several basic system parameters. The interpretation of the results is also discussed. It turns out that at the embedded epochs the behavior of the second queue is closely related to two G/M/1-type Markov processes. Utilizing the specific structure of one of these Markov processes enables to interpret these results. This gives an insight into the behavior of this tandem queuing model, which can serve as a base for developing further approximations.
Recently the tandem queueing model was also applied as a model in the area of wireless networks ( [14], [10]). Another application area is manufacturing systems and their optimal control ( [15], [13]).
From the application point of view the importance of the current work lies in the ability of modeling more general service time of the first station and the deeper understanding of the behavior of this tandem queueing model due to the interpretation of the closed-form results.
The rest of the paper is organized as follows. In section 2 we introduce and explain the model. In section 3 we deal with the structural characteristics of the model. The steady-state number of customers at the embedded epochs is investigated in section 4. In section 5 we derive the steady-state sojourn time results. It is followed by the description of the computational procedure including also numerical examples in section 6. The discussion of the results and possible further step closes the paper in section 7.
2. Model description. We consider a tandem queueing model consisting of two stations. After getting service at first station each customer moves to the second one. Each of the stations has its own queue and one server unit. Both queues have infinite buffers, i.e. all customers arriving to the system will be served at both stations. The customers arriving to the first station according to Poisson process with rate λ. The customer service times at the first station are independent and identically distributed. B, B(s), b, b (2) denote the service time r.v., its LST and its first two moments, respectively. Each customer departure from the server unit of the first station implies an immediate customer arrival to the second station. The distribution of the customer service time at second station is exponential with parameter µ. The customers during their staying at first station, i.e. during their arrival, waiting, service or departure, are called 1-customers. Similarly, the customers during their staying at second station, i.e. during their arrival, waiting, service or departure, are called 2-customers. According to this the naming of the customers moving from the first station to the second one changes from 1-customers to 2-customers.
We define D-epoch as the epoch just after the customer departure from the first station and just before the arrival of the same customer to the second station. We remark here that at the D-epoch the customer moving from the first station to the second one is neither 1-customer nor 2-customer. The different epochs and the naming of the customers are shown in Fig. 1. We also define the 1-exhaustion epoch as the D-epoch at the start of the idle period at first station. The server utilization at the first station is ρ 1 = λb. We also use the notations ρ 2 = λ µ and ρ max 2 = 1 bµ . We call this tandem queueing model also as M/G/1− > ./M/1 tandem queueing model. The model is illustrated in Fig. 2.
On this tandem queueing model we impose the following assumptions: The arrival rate and the mean customer service times at both stations are positive and finite, i.e. 0 < λ < ∞, 0 < b < ∞ and 0 < 1/µ < ∞.
A.2 The arrival process and each of the customer service times are mutually independent.
A.3 The service during the service period is work conserving at both stations.
A.4 The service during the service period is non-preemptive at both stations. A.5 The customers are served in First-In-First-Out (FIFO) order.
A. 6 The model is stable. The stability condition of this tandem queueing model is well-known and it is given as This can be explained as follows. The whole system is stable if and only if both stations are stable. The first station is an M/G/1 queue, therefore it is stable if and only if ρ 1 < 1 (see e.g. in [7]). The second station is stable if and only if the average arrival rate of the 2-customers, which is the same as the average output rate of the 1-customers, does not exceed the constant service rate at second station. However the average output rate of the 1-customers equals the arrival rate to the first station due to the stability of the first station. Thus the stability condition of the second station is given as λ < µ, which is equivalent to ρ 2 < 1.
When Y denotes a matrix, then [Y] j,l stands for its j, l-th element. Similarly [y] j denotes the j-th element of vector y. When Y(z), |z| ≤ 1 is a matrix generating function, Y denotes its value at z = 1, i.e., Y = Y(1). When y(z), |z| ≤ 1 is a vector GF, y (k) denotes its k-th (k ≥ 1) factorial moment, i.e., y (k) = d k dz k y(z)| z=1 and y denotes its value at z = 1, i.e., y = y(1). We also use notation y (k) = y (k) e. The notation eig(Y) stands for the eigenvalues of matrix Y.
3. Structural characteristics of the system. In the first part of this section we give the description of the system by means of a bivariate Markov chain embedded at D-epochs. We introduce several matrices describing the evolution of the number of 2-customers during different intervals and several partial matrix PGFs. Some of these matrices serve as a building block of the probability transition matrix of the embedded Markov chain. This transition matrix has specific structure, since it is an M/G/1-type on block level and its block matrices are G/M/1-type. In the second part of the section we characterize the evolution of the system in terms of vector GFs and matrix generating functions. This leads to the governing equation of the system for the steady-state vector GF of the number of customers at both stations.
3.1. Description of the system at D-epochs. We describe the system and its evolution at D-epochs. Let N 1 ( ) and N 2 ( ) be the number of customers at the first and at the second station at the -th D-epoch, for > 0, respectively. The sequence {(N 1 ( ), N 2 ( )), > 0} is a bivariate homogeneous embedded Markov chain on the state space ({0, 1, . . .}, {0, 1, . . .}). We say that the chain is in state (i, k) when N 1 ( ) = i and N 2 ( ) = k. We also say that the chain is in level i when N 1 ( ) = i. Let p i,k (j, n) denote the probability of transition from state (i, k) to state (j, n) in this Markov chain, i.e.
In the following we consider the structure of the transition probability matrix of the above defined Markov chain. We define γ u as the probability of u 2-customer departures during the idle period of the first station. Due to the Poisson arrival process this idle period is exponentially distributed. Using it the probability γ u can be expressed as Furthermore we define the probability γ * u as We define the ∞ × ∞ matrix Π − as Matrix Π − describes the transition in number of 2-customers during the idle period of the first station.
It follows from the definitions of γ * u and γ u that matrix Π − is stochastic, since it is nonnegative and where e is an ∞ × 1 column vector having all elements equal to one. We define β r,u as the joint probability of r 1-customer arrivals and u 2-customer departures during a 1-customer service time. The joint probability β r,u can be expressed as Moreover we also define the joint probability β * r,u as We define the ∞ × ∞ matrix Π d r as The (k, n)-th element of matrix Π d r describes the joint probabilities of r 1customer arrival and the k ⇒ n transition in number of 2-customers during the first 1-customer service time after the idle period of the first station.
Let a r denote the probability of r 1-customer arrivals during a 1-customer service time for r ≥ 0. Thus a r is given by and it can be expressed as The definition of β * r,u and (3) imply that β r,u e = a r e, r ≥ 0.
Next we define the ∞ × ∞ matrix Π (0) r as The (k, n)-th element of matrix Π (0) r describes the joint probabilities of the k ⇒ n transition in number of 2-customers from the 1-exhaustion epoch to the D-epoch at the end of the next 1-customer service and r 1-customer arrival during that 1customer service. It can be seen from their definitions that matrix Π − has G/M/1type structure, while matrix Π d r is lower triangular. It follows from (5) that Π (0) r has also G/M/1-type structure. Using (2), (4) and (5) yields We also define the ∞ × ∞ matrix Π + r as The (k, n)-th element of matrix Π + r describes the joint probabilities of the k ⇒ n transition in number of 2-customers during a 1-customer service time and r 1customer arrival during that 1-customer service. It can be seen from the structure of matrix Π + r that it has G/M/1-type structure. By using the definition of β * r,u and (3) the row sums of Π + r can be expressed as β r,u e = a r e, r ≥ 0.
We define the special Toeplitz matrix Ω 1 as , and the matrix E 0 as Observe that, using also (3), these special matrices connect the matrices Π d r and Π + r as Furthermore (5) and (6) imply that Let Π be the transition probability matrix of the above defined Markov chain in terms of block matrices, in which the (k, n)-th element of the (i, j)-th block matrix equals p i,k (j, n) for i, j, k, n ≥ 0. It follows from the interpretation of matrices Π (0) r that the (0, j)-th block matrix of Π, for j ≥ 0, equals Π (0) j . The (i, j)-th block matrix of Π, for i ≥ 1 and j ≥ 0, describes the transitions, in which j − i + 1 1-customers arrive during a 1-customer service time. Hence the (i, j)-th block matrix of Π, for i ≥ 1 and j ≥ 0, equals Π + j−i+1 . Therefore the matrix Π is given as Remark 1. The transition matrix of the Markov chain embedded at D-epochs, Π, has M/G/1-type structure on block level and its every block matrix has G/M/1-type structure.
We define the partial matrix PGF Π d (z) as

ZSOLT SAFFER AND WUYI YUE
The definition of Π d r implies that the matrix Π d (z) can be expressed as Using (8), (3) as well as several definitions gives Setting z = 1 in (8) and (9) shows that matrix Π (0) is stochastic, since it is nonnegative and we have Furthermore we define the partial matrix PGF Π The matrix Π + (z) can be expressed as

Using (3) as well as several definitions yields
Setting z = 1 in Π + (z) gives matrix Π + . It follows from the definitions of Π + r , β u (z) and β * u (z) that the matrix Π + can be expressed as Setting z = 1 in (11) shows that matrix Π + is stochastic, i.e. its every element can be interpreted as probability and

The definitions of
Finally setting z = 1 in (12) yields 3.2. The evolution of the system. We define the 1 × ∞ row vector q i , for i ≥ 0, by its k-th element as The vector q i represents the steady-state joint probabilities of the number of customers at both stations for the case when the number of 1-customers in the system is i.
We also define the corresponding steady-state vector GF q(z) as We also introduce a notation for the corresponding steady-state PGF, which is the PGF of the number of 1-customers, as Furthermore we define the 1 × ∞ row vector m as the conditional steady-state probability vector of the number of customers at the second station at 1-exhaustion epoch, i.e. given that the number of 1-customers is 0.
This definition implies that the joint probability [q 0 ] k equals [m] k P {the first station is idle at a D-epoch} for k ≥ 0. This probability equals the reciprocal of the mean number of customers served during a service period at station 1. This is 1 1−ρ1 , since the first station has an M/G/1 queue and its stochastic behavior is independent of the one of the second station. Using it we get the following relationship between q 0 and m Theorem 3.1. The governing equation of the stable M/G/1− > ./M/1 tandem queueing model satisfying assumptions A.1 -A.6 is given in terms of q(z) and m as where I denotes the ∞ × ∞ identity matrix.
Proof. Writing the equilibrium equation of the Markov chain in partitioned form yields Multiplying both sides of (16) by z i , summing over i ≥ 0 and rearrangement leads to Using the definitions of q(z) and Π + (z) in (17) gives Applying (14) in (18) and rearranging it leads to The statement of the theorem comes by applying (12) in (19) and rearranging it.
Remark 2. The structure of relation (19) is similar to that of the relation of the vector GF of the number of customers at customer departure epoch in the BMAP/G/1 queue (see e.g. in [11]). In the M/G/1− > ./M/1 tandem queueing model the number of 2-customers plays similar role as the phase of BMAP in the BMAP/G/1 queue. 4. The steady-state number of customers at D-epoch. In this section we derive several formulas starting from the governing equation of the model. The form of (15) implies that all of these formulas can be given in terms of m. Therefore first we give a method for determining the vector m. It enables the derivation of the steady-state vector GF of the number of customers at both stations, q(z), which is the main result of the paper. Then the derivation of the steady-state probability vector of the number of 2-customers at D-epochs, q, follows. In the last part of the section we provide also the expression of the vector factorial moments of the steady-state number of 1-customers at D-epochs, q (n) for n ≥ 1.

4.1.
Determination of vector m. We define matrix G, whose (k, )-th element is given as the probability that starting from state (n + 1, k) in the Markov chain defined in subsection 3.1 the first state visited in level n is (n, ), n ∈ 0, 1, 2, . . ., k, ≥ 0. Proposition 1. In the stable M/G/1− > ./M/1 tandem queueing model satisfying assumptions A.1 -A.6 the vector m is uniquely determined by the following system of linear equations: where G can be computed from the following matrix equation: Proof. The M/G/1-type structure of matrix Π ensures that matrix G can be computed by applying standard matrix-analytic method (see e.g. in [9]). This leads to Rearranging (22) results the second statement of the theorem. The sum of elements of vector m is 1, since it is a conditional probability vector, i.e. me = 1.
In order to establish a relation for m we use again a standard matrix-analytic argument. The transition from a 1-exhaustion epoch to the next 1-exhaustion epoch brings vector m to itself. Let us consider the transitions in the number of customers at the second station during the idle period and the succeeding first customer service time at the first station. If the number of 1-customer arrivals during that first customer service time is r then the above transition can be described by matrix Π (0) r . At the end of this transition there are r number of customers at the first station, i.e. the Markov chain is in level r. Getting back from level r to level 0, i.e. to idle state of the first station, requires r consecutive level decreases. Each of these level decreases can be described by matrix G. Counting for all possible values of r results in a relation for vector m as Applying (7) in (24) yields Taking into consideration (23) and using (21) in (25) gives the first statement of the theorem.

Steady-state vector GF of the number of customers.
Lemma 4.1. Let us consider the stochastic matrix Y and a Markov chain, whose transition probabilities are given by Y. If the Markov chain is irreducible and aperiodic then the matrix (I − Y + ey) is nonsingular, where y is the steady-state probability vector of that Markov chain.
Proof. If the Markov chain is irreducible and aperiodic then the probability of each state of the Markov chain after n-step transition always converge to a limit as n → ∞, i.e. the steady-state probability vector exists and it is unique. Furthermore the steady-state probability vector is a solution of the system of linear equations yY = y.
First we consider the case, when the Markov chain is positive recurrent. We apply an indirect argument. In this case y is a distribution and it is uniquely determined by the system of linear equations y(Y − I) = 0 and ye = 1, where 0 is the row vector having all elements equal zero. Let us assume that matrix (I − Y + ey) is singular, thus there exists a row vector x = 0, which fulfills the homogenous system of linear equations x(I − Y + ey) = 0. Hence x(I − Y + ey)e = 0 also holds which together with y = 0 implies that xe = 0, i.e. the row sum of x is zero. Using this it follows from the above homogenous system of linear equations that x(Y − I) = 0 is also true. This together with the uniqueness of y imply that x must be a constant multiplication of y. In this case the row sum of x must differ from zero since x = 0 and the row sum of y is not zero. However this a contradiction as it was shown that the row sum of x is zero, therefore matrix (I − Y + ey) must be nonsingular. Now consider the other case when the Markov chain is null recurrent or transient. Solving the system of linear equations yY = y by applying recursive substitution leads to y = y lim n→∞ (Y) n . If the Markov chain is null recurrent or transient then the n-step transition probabilities to each state of the chain always converge to zero as n → ∞, i.e. lim n→∞ (Y) n = 0, where 0 is the zero matrix. It follows that the only solution of the homogenous system of linear equations y(Y − I) = 0 is y = 0 and thus the matrix (Y −I) must be nonsingular. Due to y = 0 the matrices (Y −I) and (I − Y + ey) are the same, which implies the statement for this case.
Let π * (z) be the steady-state probability vector of the Markov chain whose transition matrix is stochastic matrix Π * (z). This steady-state probability vector exists and it is unique, since the above Markov chain has one irreducible class of aperiodic states. This can be seen from the structure of matrix Π * (z). The steadystate probability vector, π * (z), can be uniquely determined by π * (z)Π * (z) = π * (z) and π * (z)e = 1.
We define also the matrices Ψ(z) and Φ(z) as 1. The steady-state vector GF q(z) is given as a sum of two terms.
• The first term consists of multiplicator terms, where the first multiplicative term, q(z) is the steady-state PGF of the number of 1-customers and the second multiplicative term, π * (z) depends only on matrix GF We define the vector r(z) for the r.h.s. of (30) as Applying (31) in (30) gives q(z)(I − Π * (z)) = r(z).
Utilizing the irreducibility and aperiodicity of the Markov chain whose transition matrix is stochastic matrix Π * (z), Lemma 4.1 can be applied. Hence the matrix (I − Π * (z) + eπ * (z)) is nonsingular for real z in the region 0 < z ≤ 1. Utilizing it we multiply both sides of (33) by (I − Π * (z) + eπ * (z)) −1 Using π * (z)(I − Π * (z) + eπ * (z)) −1 = π * (z) and (13) in (34) leads to We remark here that q(z) is the PGF of the number of 1-customers in the system at D-epochs, which is the PGF of the number of customers in the M/G/1 queue due to the stochastic independency of the first station. Applying (31) in (35) and rearrangement results in This can be rewritten by using the matrices Ψ(z) and Φ(z) as Applying recursive substitution of the above expression of q(z) into the r.h.s. of (37) yields In order to investigate the convergence of the infinite matrix sum ∞ k=0 (Ψ(z)) k now we consider the eigenvalue of matrix Ψ(z) with the highest absolute value on a specific real interval. It follows from the non-singularity of matrix (I − Π * (z) + eπ * (z)) that the homogenous equation x(I−Π * (z)+eπ * (z)) = 0x has no non-trivial solution. Hence the value 0 is not an eigenvalue of matrix (I − Π * (z) + eπ * (z)). Let κ * denote the absolute value of the eigenvalue of matrix (I − Π * (z) + eπ * (z)) with the smallest absolute value. It follows that the absolute value of the eigenvalue of matrix (I − Π * (z) + eπ * (z)) −1 with the highest absolute value is 1 κ * . On the other hand B(λ − λz) goes to 1 as real-valued 0 < z ≤ 1 goes to 1. It implies that in a properly chosen (z 0 , 1) interval the expression B(λ−λz)−z B(λ−λz) can take any small value. Therefore for a given positive real number κ * there exists a real z 0 < 1 so that for every z 0 < z < 1 the relation B(λ−λz)−z B(λ−λz) < κ * holds. Taking all these into account and using the definition of Ψ(z) we conclude that there exists a real z 0 < 1 so that for every z 0 < z < 1 the absolute value of the eigenvalue of matrix Ψ(z) with the highest absolute value is less than κ * 1 κ * = 1. In other words ∃z 0 < 1 for which | eig(Ψ(z)) |< 1, z 0 < z < 1, z 0 , z ∈ R.
It follows that for z 0 < z < 1 the infinite matrix sum ∞ k=0 (Ψ(z)) k converges to a finite matrix and thus the following relations hold: Applying (39) in (38) results in the expression of the statement, (29). So far we showed the validity of the expression (29) for real z in an open real interval (z 0 , 1). Now we apply the concept of the analytic continuation in the complex analysis to extend the validity scope of (29) to |z| ≤ 1. Any PGF is an analytic complex function thus also q(z) is analytic. Moreover also the expression on the r.h.s. of (29) is an analytic complex function. It follows that they are the same in the whole convergence region of the vector GF q(z), hence the formula (29) is valid in the whole | z |≤ 1 region.

4.3.
Steady-state probability vector of the number of 2-customers. Let π + be the steady-state probability vector of the Markov chain, whose transition matrix is stochastic matrix Π + . It can be uniquely determined by π + Π + = π + and π + e = 1. 1. The steady-state probability vector q is given as a sum of two terms, in which • the first term is the steady-state probability vector of the number of customers at the customer arrival epochs of the virtual G/M/1 queue with interarrival time B and mean service time 1 µ (characterized by matrix Π + ) and • the second term represents the dependency on the difference of two G/M/1-type structure matrices, (Π (0) − Π + ). 2. The steady-state probability vector of the number of customers at the customer arrival epochs of the above virtual G/M/1 queue, π + , can be given in dependency of the stability of that virtual G/M/1 queue as • If ρ max 2 < 1 then the above virtual G/M/1 queue is stable and the probability vector π + represents a geometrical distribution, for which where σ + is the only root of the equation in the region 0 < σ + < 1.

ZSOLT SAFFER AND WUYI YUE
• If ρ max 2 ≥ 1 then the above virtual G/M/1 queue is unstable and the probability vector π + is given as 3. Furthermore the expression of q is given by Proof. The third statement of the corollary comes by setting z = 1 in (29). The steady-state probability vector of the number of customers at the customer arrival epochs of the above virtual G/M/1 queue is π + , since the transition matrix of the characterizing Markov chain is Π + . The first statement of the theorem comes from the structure of (43) and from the above consideration. The second statement is well-known from the theory of standard G/M/1 queue (see e.g. in [7]). If ρ max 2 ≥ 1, then the Markov chain whose transition matrix is stochastic matrix Π + becomes instable, i.e. it becomes either null recurrent or transient, which can be checked by applying Foster's criterion on matrix Π + .

Remark 3.
As the mean length of the idle period of station 1 tends to 0 then matrix Π (0) becomes equal to matrix Π + . As a consequence of it the second term in (43) representing the dependency on the difference of these matrices vanishes.

Corollary 2.
In the stable M/G/1− > ./M/1 tandem queueing model satisfying assumptions A.1 -A.6 the vector factorial moments of the steady-state number of 1-customers is given as a sum of two terms, from which the first term is the multiplication of • the n-th factorial moment of the number of customers in an M/G/1 queue corresponding to the first station of the model and • the steady-state probability vector of the number of customers at the customer arrival epochs of the virtual G/M/1 queue with interarrival time B and mean service time 1 µ (characterized by matrix Π + ), π + . Furthermore the expression of q (n) , for n ≥ 1, is given by where 1 (con) denotes the indicator of condition "con" and for n ≥ 1 as well as i.e. q (n) is the n-th derivative of the Pollaczek-Khinchine transform formula at z = 1.
Proof. The corollary can be proved by taking the n-th derivative of (29) and setting z = 1.

Sojourn time.
In this section we derive several results for the steady-state sojourn time of a customer in the system. The section starts by introducing several definitions of different measures related to sojourn time of the customers in the system. This is followed by the derivation of the LST of the steady-state sojourn time of the customers in the system. Based on it also the steady-state mean sojourn time of the customers in the system is derived.

5.1.
Definitions. The sojourn time of a 1-customer is defined as the time elapsed from its arrival to station 1 until its D-epoch. Let W * ,1 denote the sojourn time of a 1-customer that arrives as the -th into station 1, for ≥ 1. We define the cumulative distribution function of the steady-state sojourn time of an 1-customer, The LST of the steady-state sojourn time of an 1-customer is defined as Furthermore we define the sojourn time of a customer in the system as the time elapsed from its arrival to station 1 until its departure from station 2. Let W * denote the sojourn time of a customer that arrives as the -th into the system, for ≥ 1. We define the cumulative distribution function of the steady-state sojourn time of a customer in the system, W * (t), as The LST of the steady-state sojourn time of a customer in the system is defined as The sojourn time of a 2-customer is defined as the time elapsed from the D-epoch just before its arrival to station 2 until its departure from station 2. Let W * ,2 denote the sojourn time of a 2-customer that arrives as the -th into station 2, for ≥ 1.
Using that the customer arriving as the -th into station 1 arrives to station 2 also as the -th, W * can be given as

5.2.
LST of the sojourn time. We condition the steady-state sojourn times W * ,1 and W * ,2 on the number of 1-customers at the -th D-epoch. Thus we define the conditional cumulative distribution function of the steady-state sojourn time of a 1customer, given that the number of 1-customers at the D-epoch of that 1-customer is n as The corresponding conditional LST is defined as Furthermore we define also the steady-state joint probability W * 1 (t, n) and its LST as It follows from the definition of the conditional probability that the LST W * 1 (s|n) can be expressed as where P 1 (n) is the steady-state probability that the number of 1-customers at a D-epoch is n and it is given as Similarly we define the conditional cumulative distribution function of the steadystate sojourn time of a 2-customer, given that the number of 2-customers at the D-epoch just before the arrival of that 2-customer is n as Again the corresponding conditional LST can be defined as where W * 1 (t) can be obtained by inverse LST transform of W * 1 (s), which is given as If the number of 2-customers at D-epoch just before the arrival of the tagged 2-customer is n 2 then the sojourn time of the tagged 2-customer equals the sum of the service times of the n 2 2-customers in the second station and the service time requirements of the tagged 2-customer. Note that due to the exponential service time at the second station also the 2-customer under service needs a full service time. Thus the conditional LST of the sojourn time of the 2-customers, given that the number of 2-customers at D-epoch just before the arrival of the tagged 2-customer is n 2 , can be given as (51) Using (51) and (47) in (50) yields (52) Utilizing that the sojourn time of a 1-customer departing from station 1 is independent of the number of 1-customers arriving during it, W * 1 (s, n) can be given by Here W * 1 (t) is the inverse LST transform of W * 1 (s), which is the LST of the M/G/1 queue, since the stochastic behavior of the first station is independent of the one of the second one. This can be given by means of the Pollaczek-Khinchine transform formula (see e.g. in [7]), which results in (49).
Setting s = 0 in (53) gives P 1 (n 1 ) as The statement of the theorem comes by applying (53) and (54)  respectively. Furthermore we define the conditional steady-state mean sojourn time of a 2-customer, given that the number of 2-customers at the D-epoch just before the arrival of that 2-customer is n, as Proof. Taking the first derivative of (50) at s = 0 yields The the steady-state mean sojourn time of the 1-customers is the steady-state mean sojourn time in the standard M/G/1 queue, which is given (see e.g. [7]) as Given that the number of 2-customers at a D-epoch is n 2 ≥ 0, the conditional mean sojourn time in the second station is the sum of n 2 times the mean service time of a 2-customer and the service time of the actually arriving 2-customer. This is because the remaining service time of the 2-customer has also exponential distribution due to the exponential service time at station 2. Therefore E[W * 2 |n 2 ] can be given by The corollary comes by applying (57) and (58) Note that the steady-state total number of customers differs from the steadystate total number of customers at D-epochs, since the number of 2-customers at an arbitrary epoch and at D-epochs are different. This difference is also reflected in the last term of (59). 6. The numerical solution. In this section we deal with the computational procedure. We consider how to keep the numeric computation tractable and provide the steps of the numeric computation of the steady-state mean sojourn time of this model. We describe how we model the GI service time by Coaxian-2 distribution. In the last part of this section we present numerical examples, which demonstrates the computational procedure and show the dependency of the steady-state mean sojourn time on several system parameters.
6.1. Numerical considerations. To keep the computation of the infinite dimensional vectors and matrices tractable we apply an upper limit X on the number of 2-customers in the system, i.e. k ≤ X, and an upper limit on the number of 1-customer arrivals during a 1-customer service time, i.e. r ≤ R. This results in finite number of elements in these vectors and matrices as well as finite number of matrices Π + r and Π (0) r , which results in finite number of required operations among others in computing matrix G from (21) and vector m from (20). The proper values of X and R depend on the required precision and can be determined in an iterative manner until the difference of consecutive values of e.g. the probabilities [q] k , for every k ≤ X, becomes less than the specified error. The detailed investigation of this and other existing truncation methods and the evaluation of their impact to the computational procedure is out of scope of this work and it could be a topic of future research. 6.3. Modeling GI service time. The distribution of the GI service time is modeled by a Coaxian-2 distribution (see Fig. 3). This distribution model consists of two exponentially distributed phases with parameters µ 1 and µ 2 . After expiry of the first phase the model either takes the second phase with probability p 1 or the modeling ends with probability 1 − p. To determine the parameters of the Coaxian-2 model from the mean (b) and coefficient of variation (c b ) of the GI service time we apply the following simple two-moment fit (see [12]): The above Coaxian-2 model can realize distributions with coefficient of variation in the range c b ≥ 0.5.
We selected the Coaxian-2 distribution for the numerical examples due to its generality and simplicity. However we remark here that the results for the M/G/1− > ./M/1 tandem queueing model presented in this paper are valid for any GI service time distribution. In Fig. 4 we have plotted the steady-state mean sojourn time of the customers (E[W * ]), of the 1-customers (E[W * 1 ]) and of the 2-customers (E[W * 2 ]) as a dependency of ρ 1 besides of keeping ρ 2 fixed. This is achieved by simultaneous changing of λ and µ.  It can be seen on the figure that increasing ρ 1 leads to higher E[W * 1 ], which is clear from the first term of (55). On the other hand higher ρ 1 results in lower E[W * 2 ]. One reason for it is that higher ρ 1 is achieved by higher µ, which causes the second term in (55) (= E[W * 2 ]) directly to be less, since µ arises in its denominator. Another reason is that higher the ρ 1 lower the weight of the second term in (43), (1 − ρ 1 ), causing less probabilities in q, which also results in less E[W In Fig. 5 we have plotted the steady-state mean sojourn time of the customers (E[W * ]), of the 1-customers (E[W * 1 ]) and of the 2-customers (E[W * 2 ]) as a dependency of ρ 2 besides of keeping ρ 1 fixed. This is achieved by changing of µ. It is not surprising that changing ρ 2 has no effect on E[W * 1 ], since it depends only on ρ 1 , which is kept fixed. Thus for higher values of ρ 2 , as the figure shows, E[W * ] is dominated by E[W * 2 ]. Higher value of ρ 2 is achieved by less µ, which causes E[W * 2 ] directly to be higher, since µ arises in the denominator of this term.  Table 1 gives an overview about the behavior of the M/G/1− > ./M/1 tandem queueing model in limit ranges of the system parameters in terms of the necessary terms to compute E[W * ].

Limit
First term in (55) Second term in (55) ranges b + λb (2) 2(1−ρ1)  Table 1. Necessary terms to compute E[W * ] in different limit ranges of the system parameters The case µ 1 and ρ 1 → 1 and ρ max 2 > 1 is not considered, since it leads to an unstable system. This is because ρ 1 → 1 and ρ max 2 > 1 implies ρ 2 > 1, which follows from the following relation: Based on the behavior of the M/G/1− > ./M/1 tandem queueing model in limit ranges of the system parameters presented in Table 1 an approximation can be developed to compute E[W * ] with less computational effort. This is a potential future research task.