# A TandemBMAP/G/1→•/M/N/0 Queue with Group Occupation of Servers at the Second Station

Authors: Chesoong Kim; Alexander Dudin; Valentina Klimenok; Olga Taramin
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

We consider a two-stage tandem queue with single-server first station and multiserver second station. Customers arrive to Station 1 according to a batch Markovian arrival process (BMAP). A batch may consist of heterogeneous customers. The type of a customer is determined upon completion of a service at Station 1. The customer's type is classified based on the number of servers required to process the request of the customer at Station 2. If the required number of servers is not available, the customer may leave
the system forever or block Station 1 by waiting for the required number of servers. We determine the
stationary distribution of the system states at embedded epochs and derive the Laplace-Stieltjes transform
of the sojourn time distribution. Some key performance measures are calculated, and illustrative numerical
results are presented.


---

## Body

## 1. Introduction

Queueing networks are widely used in capacity planning and performance evaluation of computer and communication systems, service centers, and manufacturing systems among several others. Some examples of their application to real systems can be found in [1]. Tandem queues can be used for modeling real-life two-node networks as well as for the validation of general decomposition algorithms in networks (see, e.g., [2, 3]). Thus, tandem queueing systems have found much interest in the literature. An extensive survey of early papers on tandem queues can be seen in [4]. Most of these papers are devoted to exponential queueing models. Over the last two decades or so, the efforts of many investigators in tandem queues were in weakening the distribution assumptions on the service times as well as on the arrivals. In particular, the arrival process should be able to capture any correlation and burstiness that are commonly seen in the traffic of modern communication networks [3]. Such an arrival process was introduced in [5] and ever since this process is referred to as a batch Markovian arrival process (BMAP). In this paper, we deal with a tandem queue under the assumption that the customers arrive according to a BMAP.Tandem queues with theBMAP input were considered in [6–10]. The papers [6, 7, 10] are devoted to the MAP/PH/1→•/G/1 system with blocking. In [8], the tandem queues BMAP/G/1/N→•/PH/1/M-1 with losses are studied. The tandem-queue of the BMAP/G/1→•/PH/1/M-1 type with losses and feedback has been studied in [9].In the present paper, we consider aBMAP/G/1→•/M/N/0 tandem queueing system where the (possibly heterogeneous) customers arrive in batches of random sizes to Station 1. Here the customers receive service individually, and upon completion of a service the customer’s type is determined. This type identification is necessary to determine the nature of service, if any, offered at Station 2. The customer’s type is classified based on the number of servers (resources) required to process the request of the customer. The simultaneous initiation or occupation of several servers to a customer’s request is typical for the so-called nonelastic traffic in communication networks. If the required number of servers is not available at that instant of the request, the customer either leaves the system forever, or awaits until the requirement is met through the release of the sufficient number of servers. In the latter case, Station 1 will be blocked.Possible applications of the tandem queue under study lie in the modeling of the distributed server application or web server application, see, for example, [11]. Station 1 is interpreted as an authentication or an access step while Station 2 represents the computing step or data base server if the processing of a job is produced by several parallel threads. This tandem queue can model also multiaddress transmission of information. Station 1 is interpreted as a transmission channel while Station 2 regulates the transmission rate by providing necessary transmission windows (timers that are switched on at the moment of a message transmission and switched off when the receipt of this message is acknowledged or time-out expires). The performance evaluation of wireless IP networks providing heterogeneous multimedia services with different QoS demands, (see [12, Chapter 8]), is the other possible application of the model under study.In this paper, we derive the stability condition of the model under study, and briefly touch calculation of the stationary distribution of the system states at the service completion epochs at Station 1 and calculation of the system performance measures. Furthermore, we derive the Laplace-Stieltjes transform of the virtual and the actual sojourn time distributions at both stations and in the whole system. The procedures for calculation of the moments of the virtual sojourn time distribution and the mean actual sojourn time are discussed. Some numerical results illustrating the behavior of the system characteristics are presented. The problem of optimal design is numerically investigated.To the best of our knowledge, the results of our paper are novel even for the case of homogeneous customers. The most important and valuable, from the mathematical point of view, result concerns the sojourn time distribution. Previously, the sojourn time distribution in tandem queues withMAP input was considered only in [13, 14]. There, the service time distribution at both the single-server stations is of phase type which allows the authors to model the sojourn time as the time until absorption in suitably defined quasi-birth-and-death processes and continuous-time Markov chains. Because we assume general service time distribution at Station 1, we need to analyze a more complicated stochastic process.The rest of the paper is organized as follows. In Section2, the mathematical model is described. In Section 3, the results concerning the stationary distribution of the embedded Markov chain in Station 1 service completion epochs are presented. In Section 4, we focus on the analysis of the virtual and actual sojourn time distributions and their moments. In Section 5, the numerical results are presented. The paper is concluded with Section 6. Appendices contain auxiliary results, proofs, and formulas useful for computations.


## 2. The Mathematical Model

We consider a tandem queue consisting of two stations, say, Station 1 and Station 2. We assume that there is no buffer between the two stations. Station 1 is represented by theBMAP/G/1 queue. That is, the arrivals to Station 1 are described by a BMAP. The BMAP is defined by the underlying process νt,t≥0, which is an irreducible continuous time Markov chain with state space {0,…,W}, and with the matrix generating function D(z)=∑k=0∞Dkzk,|z|≤1. Arrivals occur only at epochs of the jumps in the underlying process νt,t≥0. The intensities of the transitions of the process νt accompanied by a batch of size k are defined by the matrices Dk,k≥0. The matrix D(1) is the infinitesimal generator of the process νt. The stationary distribution vector θ of this process satisfies the equations θD(1)=0,θe=1, where e is a column vector consisting of 1′s, and 0 is a row vector of 0’s.The average intensityλ (fundamental rate) of the BMAP is given by λ=θD′(z)|z=1e. We assume that λ<∞. The average intensity λb of group arrivals is defined by λb=θ(-D0)e. The coefficient of variation, cvar, of intervals between successive group arrivals is defined by cvar2=2λbθ(-D0)-1e-1. The coefficient of correlation ccor of the successive intervals between group arrivals is given by ccor=(λbθ(-D0)-1(D(1)-D0)(-D0)-1e-1)/cvar2. For more information about the BMAP and related research see, for example, [5, 15].All arriving customers enter into Station 1. The successive service times of customers at Station 1 are independent random variables with general distributionB(t), Laplace-Stieltjes transform β(s)=∫0∞e-stdB(t), and finite first moment b1=∫0∞tdB(t).After receiving a service at Station 1, the customer proceeds to Station 2. At this station, there areN identical servers. Each of these servers offers services that are exponentially distributed with parameter μ. Customers are heterogeneous with respect to the number of servers that are required to process a customer at Station 2. With probability qm,qm≥0,m=0,N¯,∑m=0Nqm=1, the customer will require exactly m servers to provide a service at Station 2 and will be called typemcustomer. Here and in the sequel, notation such as m=0,N¯, means that m assumes values from the set {0,1,…,N}. Note that customers who are all in the same batch (at the time of arriving) may belong to different types after receiving service at Station 1. Type 0 customer leaves the system for good after the service at Station 1. We assume that q0≠1. Otherwise, the queue under consideration will be reduced to the BMAP/G/1 queue which has been studied extensively.If the customer is of typem, m=1,N¯, and the required number of servers is available, the customer’s service will begin immediately. Each of these m servers processes the customer’s request independently of the others, and furthermore any server who becomes free after completing his/her share of the processing will be available to process waiting or future customers’ requests.If the required number of servers is not available, with probabilityγ,0≤γ≤1, the customer will choose to leave the system for good and with probability 1-γ will decide to wait until the required number of servers is available. In the latter case, the customer will block Station 1 since we assume that there is no buffer between the two stations. Such an assumption of blocking and loss will allow us to unify these two classes of models which are studied separately in the literature.In the following, we are interested in the steady state analysis of the model under study. For further use in the sequel, we introduce the following notation:(i)
I is an identity matrix of appropriate dimension;(ii)
⊗ and ⊕ are symbols of the Kronecker product and sum of matrices;(iii)
D̃k=IN+1⊗Dk,k≥0,D̃(z)=∑k=0∞D̃kzk,|z|≤1;(iv)
P(j,t),j≥0, is a matrix function defined by the expansion ∑j=0∞P(j,t)zj=eD(z)t;(v)
F(t)=(Fr,r′(t))r,r′=0,N¯, where Fr,r′(t)=0 for r≤r′ and for r>r′,Fr,r′(t) is the generalized Erlang distribution function with the Laplace-Stieltjes transform fr,r′(s)=∏l=r′+1rlμ(lμ+s)-1;(vi)
Qm,m=1,4¯, are square matrices: Q2=diag⁡{∑m=N-r+1Nqm,r=0,N¯},
(2.1)Q1=(q0q1…qN0q0…qN-1⋮⋮⋱⋮00…q0),Q3=(0…0qN0…0qN-1⋮⋱⋮⋮0…0q0),Q4=(q0q1…qN-1qNq0q1…qN-10⋮⋮⋱⋮⋮q00…00);(vii)
Q̃m=Qm⊗IW¯,m=1,3¯,W¯=W+1;(viii)
Q̂=Q̃1+γQ̃2+(1-γ)∫0∞(dF(t)⊗eD0t)Q̃3,Q=Q1+γQ2+(1-γ)EQ3;(ix)
(2.2)E=(00⋯0010…00⋮⋮⋱⋮⋮11⋯10),Î=(10⋯000⋯0⋮⋮⋱⋮00⋯0),Ĩ=I-Î,ê=(1,0,…,0).


## 3. The Stationary Distribution of the Embedded Markov Chain

Lettn denote the time of the nth service completion at Station 1. Consider the process ξn={in,rn,νn}, n≥1, where in,in≥0, is the number of customers at Station 1 (not counting the blocked customer, if any) at epoch tn+0; rn,rn=0,N¯, is the number of busy servers at Station 2 at epoch tn-0; νn,νn=0,W¯, is the state of the BMAP at epoch tn.It is easy to verify that the processξn={in,rn,νn}, n≥1, is a Markov chain. Enumerating the states of this Markov chain in lexicographic order, and denoting by Pl,k,l,k≥0, the square matrix of order (W+1)(N+1) governing the transition probabilities of the chain from the set of states {l,·,·} to the set {k,·,·}, the following lemma gives the entries of the transition probability matrix of the Markov chain ξn.Lemma 3.1.
The transition probability matrix of the chainξn,n≥1, has the following block structure:
(3.1)P=(Pl,k)l,k≥0=(C0C1C2C3⋯Y0Y1Y2Y3⋯0Y0Y1Y2⋯00Y0Y1⋯⋮⋮⋮⋮⋱),
where
(3.2)Ci=∑k=1i+1[-Q̂(Δ⊕D0)-1D̃k+(1-γ)FkQ̃3]Ωi-k+1,Yi=(Q̃1+γQ̃2)Ωi+(1-γ)∑k=0iFkQ̃3Ωi-k,Ωj=∫0∞eΔt⊗P(j,t)dB(t),Fj=∫0∞dF(t)⊗P(j,t),j≥0,Δ=(000⋯00μ-μ0⋯0002μ-2μ⋯00⋮⋮⋮⋱⋮⋮000⋯Nμ-Nμ).Proof.
First, we write the transition probability matricesCi,Yi in block forms as Ci=(Ci(r,r′))r,r′=0,N¯,Yi=(Yi(r,r′))r,r′=0,N¯, where the blocks Ci(r,r′), Yi(r,r′) correspond to the transitions of the number of busy servers from r to r′ at Station 2.
Denote byδr,r′(t) the probability that during the time interval of the length t the number of busy servers at Station 2 decreases from r to r′ conditioned on the fact that none arrive from Station 1.
For use in the sequel, we register the following probabilistic interpretations of the matrices.
The(ν,ν′)th entry of the matrix P(j,t) gives the probability that j customers arrive in the BMAP during the interval (0,t] and the state of the BMAP at epoch t is ν′ given ν0=ν.
The(ν,ν′)th entry of the matrix ∫0∞δr,r′(t)P(j,t)dB(t) gives the probability that during the service time of a customer at Station 1, exactly j customers arrive, the number of busy servers at Station 2 decreases from r to r′ and the BMAP has changed from ν to ν′.
The(ν,ν′)th entry of the matrix ∫0∞δr,r′(t)eD0tDkdt gives the probability that at an arbitrary time instant with the number of busy servers at Station 2 equal to r, and the BMAP in state ν, the first batch to arrive is of size k; soon after that instant, the number of busy servers at Station 2 equals r′ and the BMAP is in state ν′.
Fr,r′(t) is the distribution function of the time interval during which the number of busy servers at Station 2 decreases from r to r′ conditioned on the fact that none arrive to this station. Then the (ν,ν′)th entry of the matrix ∫0∞P(j,t)dFr,r′(t) defines the probability that exactly j customers arrive with the BMAP moving from ν to ν′ and that the number of busy servers at Station 2 decreases from r to r′ during that time interval.
From the above probabilistic interpretations, analyzing the one-step transitions of the chainξn with a careful analysis of the service completion epochs, in which the customers may get lost due to lack of servers or wait (and thus block Station 1) until enough servers are available, we obtain the following expressions for the matrices Ci(r,r′),Yi(r,r′),i≥0:
(3.3)Ci(r,r′)=∑m=0N-rqm∑l=r′r+m∫0∞δr+m,l(t)eD0tdt∑k=1i+1Dk∫0∞δl,r′(t)P(i-k+1,t)dB(t)+γ∑m=N-r+1Nqm∑l=r′r∫0∞δr,l(t)eD0tdt∑k=1i+1Dk∫0∞δl,r′(t)P(i-k+1,t)dB(t)+(1-γ)∑m=N-r+1Nqm[∫0∞eD0tdFr,N-m(t)∑l=r′N∫0∞δN,l(t)eD0tdt×∑k=1i+1Dk∫0∞δl,r′(t)P(i-k+1,t)dB(t)+∑k=1i+1∫0∞P(k,t)dFr,N-m(t)×∫0∞δN,r′(t)P(i-k+1,t)dB(t)],Yi(r,r′)=∑m=0N-rqm∫0∞P(i,t)δr+m,r′(t)dB(t)+∑m=N-r+1Nqm[γ∫0∞δr,r′(t)P(i,t)dB(t)+(1-γ)∑k=0i∫0∞P(k,t)dFr,N-m(t)∫0∞δN,r′(t)P(i-k,t)dB(t)].
In order to arrive at equations (3.2) from (3.3), we use the matrix notations introduced above and the relation: (δr,r′(t))r,r′=0,N¯=eΔt, which follows from the fact that under the case when none arrive to Station 2, the process rt governing the number of busy servers at this station is Markovian with generator Δ.It is easy to see that the Markov chainξn belongs to the class of M/G/1 type Markov chains, see [16]. We can use this fact to derive the ergodicity condition and calculate the stationary distribution of the chain.LetC(z)=∑i=0∞Cizi,Y(z)=∑i=0∞Yizi,|z|≤1, be the generating functions of the transition probability matrices Ci and Yi,i≥0.Corollary 3.2.
The matrix generating functionsC(z),Y(z) can be written as
(3.4)C(z)=1z[-Q̂(Δ⊕D0)-1(D̃(z)-D̃0)+(1-γ)(F(z)-F0)Q̃3]Ω(z),(3.5)Y(z)=[Q̃1+γQ̃2+(1-γ)F(z)Q̃3]Ω(z),
where
(3.6)Ω(z)=∑n=0∞Ωnzn=∫0∞eΔt⊗eD(z)tdB(t),F(z)=∑n=0∞Fnzn=∫0∞dF(t)⊗eD(z)t.Theorem 3.3.
The necessary and sufficient condition for ergodicity of the Markov chainξn,n≥1, is the fulfillment of the inequality
(3.7)ρ=λ[b1+(1-γ)∑r=1Nϑr∑m=N-r+1Nqm∑l=N-m+1r(lμ)-1]<1.
Here (ϑ1,…,ϑN) is a part of the vector ϑ=(ϑ0,…,ϑN), which is the unique solution to the system
(3.8)ϑQB*(0)=ϑ,ϑe=1,
where B*(s)=∫0∞e-steΔtdB(t).Proof.
It can be verified that the matrixY(1) is irreducible. Hence, from [16], the necessary and sufficient condition for ergodicity of the chain ξn is the fulfillment of the inequality
(3.9)xY′(1)e<1,
where the vector x is the unique solution of the system
(3.10)xY(1)=x,xe=1.
The theorem will be proven if we show that inequality (3.9) is equivalent to inequality (3.7).
Let the vectorx be of the form (3.11)x=ϑ⊗θ.
By the direct substitution into the system (3.10), where Y(1) is calculated using (3.5), we verify that such a vector provides the unique solution of this system. Differentiating (3.5) at the point z=1 and substituting the resulting expression for Y′(1) and the vector x of form (3.11) into the inequality (3.9), we get
(3.12)ρ=λ[b1+(1-γ)ϑ∫0∞tdF(t)Q3e]<1.
The stated expression in (3.7) follows from (3.12) and the expression for ∫0∞tdF(t) given in (C.3)-(C.4) (see Appendix C).Remark 3.4.
The inequality (3.7) is intuitively clear on noting that the vector ϑ gives the stationary distribution of the number of busy servers at Station 2 at the service completion epochs at Station 1 given the latter station works non-stop. Then (1-γ)∑r=1Nϑr∑m=N-r+1Nqm∑l=N-m+1r(lμ)-1 defines the average blocking time of Station 1 under overload condition and ρ is the system load.Remark 3.5.
In a majority of queueing systems, the system load depends only on the first moment of the service time distribution. In the model under study, the value ofρ depends not only on the first moment b1 of the service time distribution at Station 1, but also on the shape of this distribution. In particular, ρ depends on variance of the service time. This fact is illustrated in Table 1 in Section 5.Table 1
The value of the system loadρ for different value of the mean service time and service time variation.
cvar=0cvar=1cvar=5cvar=9.95b1=0.10.427540.430160.455570.47940b1=0.20.474600.482280.536770.57784b1=0.30.531730.545150.622430.67725b1=0.40.596980.615900.710100.77695b1=0.50.668790.692470.799020.87678b1=0.60.745790.773390.888830.97667b1=0.70.826870.857630.979371.07663b1=0.80.911160.944401.070511.17655In what follows, we assume that the inequality (3.7) holds true.Denote the stationary state probabilities of the Markov chainξn={in,rn,νn} by π(i,r,ν),i≥0,r=0,N¯,ν=0,W¯. Introduce the notation for the row vectors of these probabilities(3.13)π(i,r)=(π(i,r,0),π(i,r,1),…,π(i,r,W)),πi=(π(i,0),π(i,1),…,π(i,N)),i≥0.
Let also Π(z)=∑i=0∞πizi,|z|≤1, be the vector generating function of vectors πi,i≥0. To compute these vectors as well as the vectors Π(1) and Π′(1), known algorithms, see, for example, [16], can be applied.Once the stationary distribution has been computed, we can calculate some key performance measures of the system as follows.(i)
The mean number of customers at Station 1 at the service completion epochsL=Π′(1)e.(ii)
The vector of the stationary distribution of the number of busy servers at Station 2 at the service completion epoch at Station 1:r=Π(1)(IN+1⊗eW¯).(iii)
The mean number of busy servers at Station 2 at the service completion epoch at Station 1(3.14)Nbusy=rdiag⁡{r,r=0,N¯}e.(iv)
The probability that an arbitrary customer leaves the system or causes the blocking of the server at Station 1(3.15)Ploss=γΠ(1)Q̃2e,Pblock=(1-γ)Π(1)Q̃2e.(v)
The probability that the server of Station 1 is idle at an arbitrary timepidle=τ-1π0Q̂(-D̃0)-1e, where τ is the mean interdeparture time at Station 1,
(3.16)τ=b1+π0Q̂(-D̃0)-1e-(1-γ)Π(1)(IN+1⊗e)F(1)Q3e,
the matrix F(1) is defined by formula (C.3) below.(vi)
The probability that the server of Station 1 processes a customer at an arbitrary timepserve=τ-1b1.(vii)
The probability that the server of Station 1 is blocked at an arbitrary timepblock=1-pidle-pserve.


## 4. Stationary Distribution of the Sojourn Time

### 4.1. The Virtual Sojourn Time

The virtual sojourn time in the system consists of the virtual sojourn time at Station 1 and the sojourn time at Station 2. We assume that customers are served according toFIFO (first-in-first-out) discipline.For use in the sequel, we define thegeneralized service time of an arbitrary customer as the service time of this customer by the first server and the possible blocking time of the server by the previous customer.


#### 4.1.1. The Virtual Sojourn Time at Station 1

The virtual sojourn time at Station 1 consists of (i) theresidual time from an arbitrary time instant (associated with virtual customer arrival) to the next service completion epoch at Station 1 (ii) the generalized service times of customers staying in the queue at an arbitrary time, and (iii) the generalized service time of the virtual customer.First, we study theresidual time. To this end, we consider the process χt={it,mt,rt,νt,ṽt},t≥0, whose components are defined as follows: it is the number of customers in Station 1 (including the blocked customer, if any), mt takes values 0, 1, 2, respectively, based on the server at Station 1 is idle, busy, or blocked at time t, νt is the state of the BMAP, rt is the number of busy servers at Station 2 just before the service completion epoch following the time t, ṽt is the residual time from t to that service completion epoch.Using the definition of semiregenerative processes given in [17], it can be verified that the process χt is a semi-regenerative one with the embedded Markov renewal process {ξn,tn},n≥1. Let(4.1)Ṽ(i,m,r,ν,x)=lim⁡t→∞P{it=i,mt=m,rt=r,νt=ν,ṽt<x},i≥0,m=0,2¯,r=0,N¯,ν=0,W¯,x≥0,
be the stationary distribution of the process χt,t≥0.From [17], the limits in (4.1) exist if the process {ξn,tn},n≥1, is irreducible aperiodic recurrent and the value τ of the mean inter-departure time at Station 1 (given by (3.16)) is finite. All these conditions hold if inequality (3.7) is satisfied.LetṼ(i,m,x) be the row vector of the steady state probabilities Ṽ(i,m,r,ν,x) arranged according to the lexicographic order of the components (r,ν), and let ṽ(i,m,s) be the corresponding vector of the Laplace-Stieltjes transforms, that is, let ṽ(i,m,s)=∫0∞e-sxdṼ(i,m,x),i≥0,m=0,2¯.Lemma 4.1.
The vector Laplace-Stieltjes transformsṽ(i,m,s) are calculated by
(4.2)ṽ(0,1,s)=0,ṽ(i,0,s)=0,i>0,ṽ(0,2,s)=0,(4.3)ṽ(0,0,s)=-τ-1π0Q̂(Δ⊕D0)-1[B*(s)⊗IW¯],(4.4)ṽ(i,1,s)=τ-1{π0∑k=1i[-Q̂(Δ⊕D0)-1D̃k+(1-γ)FkQ̃3]×∫0∞(eΔu⊗IW¯)∫0uIN+1⊗P(i-k,y)e-s(u-y)dydB(u)+∑j=1iπj[(Q̃1+γQ̃2)∫0∞(eΔu⊗IW¯)∫0uIN+1⊗P(i-j,y)e-s(u-y)dydB(u)+(1-γ)∑k=0i-jFkQ̃3∫0∞eΔu⊗P(i-k-j,y)e-s(u-y)dydB(u)]},(4.5)ṽ(i,2,s)=τ-1(1-γ)∑j=0i-1πj(∫0∞dF(u)⊗IW¯)Q̃3×∫0ue-s(u-y)[IN+1⊗P(i-j-1,y)]dy[B*(s)⊗IW¯],i>0.Proof.
Letκj(i,m,x,t)(r,ν;r′,ν′) denote the conditional probability that, given time 0 is an instant of the service completion at Station 1 and the embedded Markov chain ξn is in the state (j,r,ν) at that time, the next service completion epoch at Station 1 occurs later than t, the discrete components of the process χt take values (i,m,r′,ν′) at time t and the continuous-time component ṽt<x.
Let us arrange the probabilitiesκj(i,m,x,t)(r,ν;r′,ν′), for fixed values i,j,m, according to the lexicographic order of the states (r,ν;r′,ν′) and form the square matrices
(4.6)K̃j(i,m,x,t)=(κj(i,m,x,t)(r,ν;r′,ν′))ν,ν′=0,W¯;r,r′=0,N¯.
Then, using the ergodic theorem for semi-regenerative processes, (see [17, Theorem  6.12]), the probability vectors Ṽ(i,m,x) can be related to the stationary distribution πj,j≥0, of the embedded Markov chain ξn,n≥1, by
(4.7)Ṽ(i,m,x)=τ-1∑j=0∞πj∫0∞K̃j(i,m,x,t)dt,i≥0,m=0,2¯.
The corresponding vector Laplace-Stieltjes transforms ṽ(i,m,s) are defined by
(4.8)ṽ(i,m,s)=τ-1∑j=0∞πj∫0∞K̃j*(i,m,s,t)dt,i≥0,m=0,2¯,
where K̃j*(i,m,s,t)=∫0∞e-sxdK̃j(i,m,x,t).
From (4.8), formulas (4.2) follow immediately when we note that K̃j*(i,m,s,t)=0 for the range of arguments {j≥0,i=0,m=1,2} and {j≥0,i>0,m=0}.
Letm=1. The lengthy but straightforward expressions for the matrices K̃j*(i,1,s,t), i>0, are presented in Appendix A. Substituting these expressions into (4.8) and after routine algebraic manipulations including rearranging the order of integration, we get formula (4.4) for the vectors ṽ(i,1,s),i>0. Similar calculations yield (4.3) and (4.5).Further, we study thegeneralized service time distribution at Station 1.LetB̂(x) be the matrix distribution function of generalized service time. More specifically, let B̂(x)=(B̂(x)r,r′)r,r′=0,N¯, where B̂(x)r,r′=P{tn+1-tn<x,rn+1=r′∣rn=r,in≠0}. Denote ℬ(s)=∫0∞e-stdB̂(t).Lemma 4.2.
The matrix Laplace-Stieltjes transform of the generalized service time distribution at Station 1 is calculated as(4.9)B(s)=[Q1+γQ2+(1-γ)F*(s)Q3]B*(s),
where F*(s)=∫0∞e-stdF(t).Proof.
To prove, we need to analyze the structure of the generalized service time. The generalized service time of a tagged customer is just the service time of the customer at Station 1 if the previous customer did not block the server of this station. In this case, the matrix Laplace-Stieltjes transform of the generalized service time distribution is calculated by(Q1+γQ2)B*(s). However, when blocking occurs, the generalized service time consists of the time during which the server is blocked by the previous customer and the service time of the tagged customer. The corresponding Laplace-Stieltjes transform is defined by (1-γ)∫0∞e-stdF(t)Q3B*(s). The stated result (4.9) follows immediately.Now we are ready to derive the equation for the vector Laplace-Stieltjes transformv1(s) of the distribution of the virtual sojourn time at Station 1. Let v1(r,ν,x) be the probability that, at an arbitrary epoch, the BMAP is in state ν, the virtual sojourn time at Station 1 is less than x, and the number of busy servers at Station 2 just before the end of the virtual sojourn time is r. Then v1(s) is defined as a vector of Laplace-Stieltjes transforms v1(r,ν,s)=∫0∞e-sxdv1(r,ν,x) written in lexicographic order.Theorem 4.3.
The vector Laplace-Stieltjes transformv1(s) satisfies the equation
(4.10)v1(s)A(s)=π0Φ(s),
where
(4.11)A(s)=sI+∑r=0∞Br(s)⊗Dr,Φ(s)=τ-1Q̂(Δ⊕D0)-1(Δ⊗IW¯-sI)[B*(s)⊗IW¯].Proof.
As mentioned above, the virtual sojourn time at Station 1 consists of the residual time from an arbitrary timet to the next service completion epoch, the generalized service times of customers that await for a service at time t, and the generalized service time of the virtual customer.
Taking into account the structure of the virtual sojourn time and using the law of total probability, we express the vector Laplace-Stieltjes transformv1(s) as follows:(4.12)v1(s)=ṽ(0,0,s)+∑i=1∞ṽ(i,1,s)[Bi(s)⊗IW¯]+∑i=1∞ṽ(i,2,s)[Bi-1(s)⊗IW¯].
Further, we multiply (4.12) by the matrix sI+∑r=0∞ℬr(s)⊗Dr, and, after some laborious calculations, get(4.13)v1(s)(sI+∑r=0∞Br(s)⊗Dr)=τ-1{π0∑i=0∞[Ci+∑j=1i+1πjYi-j+1][Bi+1(s)⊗IW¯]+π0Q̂[-(Δ⊕D0)-1(sI+D̃0)+I][B*(s)⊗IW¯]-∑j=0∞πj[Bj+1(s)⊗IW¯]}.
Multiplying the balance equations for stationary probability vectors πi of the form
(4.14)πi=π0Ci+∑l=1i+1πlYi-l+1,i≥0,
by ℬi+1(s)⊗IW¯ and summing over i we obtain
(4.15)∑i=0∞πi[Bi+1(s)⊗IW¯]=π0∑i=0∞Ci[Bi+1(s)⊗IW¯]+∑i=0∞∑j=1i+1πjYi-j+1[Bi+1(s)⊗IW¯].
Using (4.15) to simplify equation (4.13), we obtain (4.10).


#### 4.1.2. The Sojourn Time at Station 2

Letv2(s) be the column vector of the Laplace-Stieltjes transforms of the conditional sojourn time distributions at Station 2. The rth entry of this vector is the Laplace-Stieltjes transform of the sojourn time distribution of a customer at Station 2 given that the number of busy servers is equal to r just before the end of the sojourn time of this customer at Station 1.Lemma 4.4.
The vector Laplace-Stieltjes transform of the sojourn time distribution at Station 2 is given by(4.16)v2(s)=[Q4(F*(s)+I)Î+γQ2+(1-γ)F*(s)diag⁡{fr,0(s),r=N,N-1,…,0}Q3]e.Proof.
The sojourn time of a customer who requiresm servers at Station 2 consists of:(i)
the service time of the customer when at leastm servers are available at the time of the request;(ii)
zero time, if the required number of servers is not available and the customer leaves the system;(iii)
the blocking time and the service time of a customer when the required number of servers is not available and the customer awaits the release of sufficiently many servers.
Note that we assume that the service of typem is performed by m servers independently of each other and finishes when all m servers complete the service. The distribution of this service time is defined by the Laplace-Stieltjes transform fm,0(s),m=1,N¯.
Taking into account this fact together with (i)–(iii) and using the matrix notation, we obtain expression (4.16) for the vector v2(s). In the expression, the first summand corresponds to the case (i), the second and the third summands give the Laplace-Stieltjes transform under study in the cases (ii) and (iii), respectively.


#### 4.1.3. The Virtual Sojourn Time in the System

Theorem 4.5.
The Laplace-Stieltjes transform of the virtual sojourn time distribution in the system is given by(4.17)v(s)=v1(s)(IN+1⊗eW¯)v2(s),
where the vectors v1(s) and v2(s), respectively, are as given in (4.10) and (4.16).Proof.
Formula (4.17) readily follows from the structure of the virtual sojourn time in the system which consists of the virtual sojourn time at Station 1 and the sojourn time at Station 2.


### 4.2. The Actual Sojourn Time

Letv1(a)(s) and v(a)(s) be the Laplace-Stieltjes transforms of the distribution of the actual sojourn time at Station 1 and in the whole system.Theorem 4.6.
The Laplace-Stieltjes transform of the actual sojourn time distribution at Station 1 is calculated as follows:(4.18)v1(a)(s)=λ-1v1(s)∑k=0∞[Bk(s)(B(s)-I)-1⊗Dk]e.Proof.
The actual sojourn time at Station 1 of an arbitrary-tagged customer, who arrived in a group of sizek and placed at the jth position within the group, consists of (a) the actual sojourn time at Station 1 of the first customer in the group, which coincides with the virtual sojourn time at Station 1; (b) the generalized service times at Station 1 of the j-2 customers of the group who arrived with the tagged customer; (c) the generalized service time of the tagged customer at Station 1.
The vector Laplace-Stieltjes transform of the sojourn time distribution at Station 1 of the first customer of thek-size group that contains the tagged customer is evidently given by the vector v1(s)(IN+1⊗kDke/λ).
Assuming that an arbitrary customer arriving in a group of sizek is placed on the jth position with probability 1/k and using the law of total probability, we immediately obtain the following expression:
(4.19)v1(a)(s)=∑k=1∞v1(s)(IN+1⊗kDkeλ)∑j=1k1kBj-1(s)e.
After some algebraic manipulations (4.19) is reduced to (4.18).Corollary 4.7.
The Laplace-Stieltjes transform of the actual sojourn time distribution in the whole system is calculated as follows:(4.20)v(a)(s)=λ-1v1(s)∑k=0∞[Bk(s)(B(s)-I)-1⊗Dk](IN+1⊗eW¯)v2(s).


### 4.3. Moments of the Sojourn Time Distribution

The formulas for the moments of the virtual sojourn time distribution can be obtained by differentiating the expression in (4.10) at the point s=0. This requires the calculation of the derivatives of v1(s) at the point s=0. However, the matrix A(s) in (4.10) is singular at the point s=0 and calculation of the derivatives at this point is the nontrivial task. The results given below allow one to develop a procedure for calculating the required derivatives.We will use notationv1(m)(s) for the mth derivative of the vector v1(s),m≥1,and set v1(0)(s)=v1(s). Similar notations will be used for other functions of s.Theorem 4.8.
Let∫0∞tmdB(t)<∞,m=1,M+1¯, where M is an arbitrary positive integer. Then the vectors v1(m)(0),m=1,M¯, are computed recursively by
(4.21)v1(m)(0)=[(π0Φ(m)(0)-∑l=0m-1(ml)v1(l)(0)A(m-l)(0))Ĩ+1m+1(π0Φ(m+1)(0)-∑l=0m-1(m+1l)v1(l)(0)A(m+1-l)(0))eê]Ã-1,
with the initial condition
(4.22)v1(0)(0)=v1(0)=[τ-1π0Q̂(Δ⊕D0)-1(ΔB*(0)⊗IW¯)Ĩ+pidleê]Ã-1,
where Ã=A(0)Ĩ+A′(0)eê.Proof of the theorem is presented in AppendixB. In what follows we assume that ∫0∞tkdB(t)<∞,k=1,2¯.Corollary 4.9.
The mean virtual sojourn time at Station 1 is given by(4.23)v¯1={[τ-1π0Q̂(Δ⊕D0)-1((B*(0)-ΔB*′(0))⊗IW¯)+v1(0)A′(0)]Ĩ+[pidleb1+12v1(0)A′′(0)e]ê}Ã-1e,
where the vector v1(0) is defined by formula (4.22).Proof follows from the relationv¯1=-v1′(0)e and formula (4.21).Theorem 4.10.
The mean virtual sojourn time in the system is given by(4.24)v¯=v¯1+v1(0)(IN+1⊗eW¯)v¯2,
where v1(0) and v¯1 are given in (4.22) and (4.23), respectively, and v¯2 is a vector of conditional means of the sojourn time at Station 2,
(4.25)v¯2=-{Q4F(1)êT+(1-γ)[F(1)-Ediag⁡{∑l=1r(lμ)-1,r=N,N-1,…,0}]Q3e},
where F(1) is given by formula (C.3) below.Proof.
To calculate the valuev¯ we differentiate the expression in (4.17). Setting s=0 and replacing the sign, we have that
(4.26)v¯=-v1′(0)(IN+1⊗eW¯)v2(0)-v1(0)(IN+1⊗eW¯)v2′(0).
Putting s=0 in (4.16) we get v2(0)=[Q4(E+I)Î+γQ2+(1-γ)EQ3]e=e. This implies that the first term in the right-hand side of (4.26) is equal to -v1′(0)e=v¯1. Using the relation v¯2=-v2′(0) and differentiating (4.16) at the point s=0, we readily verify that v¯2 has form (4.25). This completes the proof.Theorem 4.11.
The mean actual sojourn time at Station 1 is given by(4.27)v¯1(a)=-λ-1{v1′(0)(e⊗∑k=1∞kDke)+v1(0)∑k=1∞(IN+1⊗Dke)[∑n=1k-1∑l=0n-1Bl(0)B′(0)e]}.Proof of the theorem follows from the relationv¯1(a)=-dv1(a)(s)/ds|s=0 and formula (4.19).Corollary 4.12.
The mean actual sojourn time in the system is given by(4.28)v¯a=v¯1(a)+λ-1v1(0)∑k=1∞(IN+1⊗Dke)∑n=0k-1Bn(0)v¯2.


## 4.1. The Virtual Sojourn Time

The virtual sojourn time in the system consists of the virtual sojourn time at Station 1 and the sojourn time at Station 2. We assume that customers are served according toFIFO (first-in-first-out) discipline.For use in the sequel, we define thegeneralized service time of an arbitrary customer as the service time of this customer by the first server and the possible blocking time of the server by the previous customer.


### 4.1.1. The Virtual Sojourn Time at Station 1

The virtual sojourn time at Station 1 consists of (i) theresidual time from an arbitrary time instant (associated with virtual customer arrival) to the next service completion epoch at Station 1 (ii) the generalized service times of customers staying in the queue at an arbitrary time, and (iii) the generalized service time of the virtual customer.First, we study theresidual time. To this end, we consider the process χt={it,mt,rt,νt,ṽt},t≥0, whose components are defined as follows: it is the number of customers in Station 1 (including the blocked customer, if any), mt takes values 0, 1, 2, respectively, based on the server at Station 1 is idle, busy, or blocked at time t, νt is the state of the BMAP, rt is the number of busy servers at Station 2 just before the service completion epoch following the time t, ṽt is the residual time from t to that service completion epoch.Using the definition of semiregenerative processes given in [17], it can be verified that the process χt is a semi-regenerative one with the embedded Markov renewal process {ξn,tn},n≥1. Let(4.1)Ṽ(i,m,r,ν,x)=lim⁡t→∞P{it=i,mt=m,rt=r,νt=ν,ṽt<x},i≥0,m=0,2¯,r=0,N¯,ν=0,W¯,x≥0,
be the stationary distribution of the process χt,t≥0.From [17], the limits in (4.1) exist if the process {ξn,tn},n≥1, is irreducible aperiodic recurrent and the value τ of the mean inter-departure time at Station 1 (given by (3.16)) is finite. All these conditions hold if inequality (3.7) is satisfied.LetṼ(i,m,x) be the row vector of the steady state probabilities Ṽ(i,m,r,ν,x) arranged according to the lexicographic order of the components (r,ν), and let ṽ(i,m,s) be the corresponding vector of the Laplace-Stieltjes transforms, that is, let ṽ(i,m,s)=∫0∞e-sxdṼ(i,m,x),i≥0,m=0,2¯.Lemma 4.1.
The vector Laplace-Stieltjes transformsṽ(i,m,s) are calculated by
(4.2)ṽ(0,1,s)=0,ṽ(i,0,s)=0,i>0,ṽ(0,2,s)=0,(4.3)ṽ(0,0,s)=-τ-1π0Q̂(Δ⊕D0)-1[B*(s)⊗IW¯],(4.4)ṽ(i,1,s)=τ-1{π0∑k=1i[-Q̂(Δ⊕D0)-1D̃k+(1-γ)FkQ̃3]×∫0∞(eΔu⊗IW¯)∫0uIN+1⊗P(i-k,y)e-s(u-y)dydB(u)+∑j=1iπj[(Q̃1+γQ̃2)∫0∞(eΔu⊗IW¯)∫0uIN+1⊗P(i-j,y)e-s(u-y)dydB(u)+(1-γ)∑k=0i-jFkQ̃3∫0∞eΔu⊗P(i-k-j,y)e-s(u-y)dydB(u)]},(4.5)ṽ(i,2,s)=τ-1(1-γ)∑j=0i-1πj(∫0∞dF(u)⊗IW¯)Q̃3×∫0ue-s(u-y)[IN+1⊗P(i-j-1,y)]dy[B*(s)⊗IW¯],i>0.Proof.
Letκj(i,m,x,t)(r,ν;r′,ν′) denote the conditional probability that, given time 0 is an instant of the service completion at Station 1 and the embedded Markov chain ξn is in the state (j,r,ν) at that time, the next service completion epoch at Station 1 occurs later than t, the discrete components of the process χt take values (i,m,r′,ν′) at time t and the continuous-time component ṽt<x.
Let us arrange the probabilitiesκj(i,m,x,t)(r,ν;r′,ν′), for fixed values i,j,m, according to the lexicographic order of the states (r,ν;r′,ν′) and form the square matrices
(4.6)K̃j(i,m,x,t)=(κj(i,m,x,t)(r,ν;r′,ν′))ν,ν′=0,W¯;r,r′=0,N¯.
Then, using the ergodic theorem for semi-regenerative processes, (see [17, Theorem  6.12]), the probability vectors Ṽ(i,m,x) can be related to the stationary distribution πj,j≥0, of the embedded Markov chain ξn,n≥1, by
(4.7)Ṽ(i,m,x)=τ-1∑j=0∞πj∫0∞K̃j(i,m,x,t)dt,i≥0,m=0,2¯.
The corresponding vector Laplace-Stieltjes transforms ṽ(i,m,s) are defined by
(4.8)ṽ(i,m,s)=τ-1∑j=0∞πj∫0∞K̃j*(i,m,s,t)dt,i≥0,m=0,2¯,
where K̃j*(i,m,s,t)=∫0∞e-sxdK̃j(i,m,x,t).
From (4.8), formulas (4.2) follow immediately when we note that K̃j*(i,m,s,t)=0 for the range of arguments {j≥0,i=0,m=1,2} and {j≥0,i>0,m=0}.
Letm=1. The lengthy but straightforward expressions for the matrices K̃j*(i,1,s,t), i>0, are presented in Appendix A. Substituting these expressions into (4.8) and after routine algebraic manipulations including rearranging the order of integration, we get formula (4.4) for the vectors ṽ(i,1,s),i>0. Similar calculations yield (4.3) and (4.5).Further, we study thegeneralized service time distribution at Station 1.LetB̂(x) be the matrix distribution function of generalized service time. More specifically, let B̂(x)=(B̂(x)r,r′)r,r′=0,N¯, where B̂(x)r,r′=P{tn+1-tn<x,rn+1=r′∣rn=r,in≠0}. Denote ℬ(s)=∫0∞e-stdB̂(t).Lemma 4.2.
The matrix Laplace-Stieltjes transform of the generalized service time distribution at Station 1 is calculated as(4.9)B(s)=[Q1+γQ2+(1-γ)F*(s)Q3]B*(s),
where F*(s)=∫0∞e-stdF(t).Proof.
To prove, we need to analyze the structure of the generalized service time. The generalized service time of a tagged customer is just the service time of the customer at Station 1 if the previous customer did not block the server of this station. In this case, the matrix Laplace-Stieltjes transform of the generalized service time distribution is calculated by(Q1+γQ2)B*(s). However, when blocking occurs, the generalized service time consists of the time during which the server is blocked by the previous customer and the service time of the tagged customer. The corresponding Laplace-Stieltjes transform is defined by (1-γ)∫0∞e-stdF(t)Q3B*(s). The stated result (4.9) follows immediately.Now we are ready to derive the equation for the vector Laplace-Stieltjes transformv1(s) of the distribution of the virtual sojourn time at Station 1. Let v1(r,ν,x) be the probability that, at an arbitrary epoch, the BMAP is in state ν, the virtual sojourn time at Station 1 is less than x, and the number of busy servers at Station 2 just before the end of the virtual sojourn time is r. Then v1(s) is defined as a vector of Laplace-Stieltjes transforms v1(r,ν,s)=∫0∞e-sxdv1(r,ν,x) written in lexicographic order.Theorem 4.3.
The vector Laplace-Stieltjes transformv1(s) satisfies the equation
(4.10)v1(s)A(s)=π0Φ(s),
where
(4.11)A(s)=sI+∑r=0∞Br(s)⊗Dr,Φ(s)=τ-1Q̂(Δ⊕D0)-1(Δ⊗IW¯-sI)[B*(s)⊗IW¯].Proof.
As mentioned above, the virtual sojourn time at Station 1 consists of the residual time from an arbitrary timet to the next service completion epoch, the generalized service times of customers that await for a service at time t, and the generalized service time of the virtual customer.
Taking into account the structure of the virtual sojourn time and using the law of total probability, we express the vector Laplace-Stieltjes transformv1(s) as follows:(4.12)v1(s)=ṽ(0,0,s)+∑i=1∞ṽ(i,1,s)[Bi(s)⊗IW¯]+∑i=1∞ṽ(i,2,s)[Bi-1(s)⊗IW¯].
Further, we multiply (4.12) by the matrix sI+∑r=0∞ℬr(s)⊗Dr, and, after some laborious calculations, get(4.13)v1(s)(sI+∑r=0∞Br(s)⊗Dr)=τ-1{π0∑i=0∞[Ci+∑j=1i+1πjYi-j+1][Bi+1(s)⊗IW¯]+π0Q̂[-(Δ⊕D0)-1(sI+D̃0)+I][B*(s)⊗IW¯]-∑j=0∞πj[Bj+1(s)⊗IW¯]}.
Multiplying the balance equations for stationary probability vectors πi of the form
(4.14)πi=π0Ci+∑l=1i+1πlYi-l+1,i≥0,
by ℬi+1(s)⊗IW¯ and summing over i we obtain
(4.15)∑i=0∞πi[Bi+1(s)⊗IW¯]=π0∑i=0∞Ci[Bi+1(s)⊗IW¯]+∑i=0∞∑j=1i+1πjYi-j+1[Bi+1(s)⊗IW¯].
Using (4.15) to simplify equation (4.13), we obtain (4.10).


### 4.1.2. The Sojourn Time at Station 2

Letv2(s) be the column vector of the Laplace-Stieltjes transforms of the conditional sojourn time distributions at Station 2. The rth entry of this vector is the Laplace-Stieltjes transform of the sojourn time distribution of a customer at Station 2 given that the number of busy servers is equal to r just before the end of the sojourn time of this customer at Station 1.Lemma 4.4.
The vector Laplace-Stieltjes transform of the sojourn time distribution at Station 2 is given by(4.16)v2(s)=[Q4(F*(s)+I)Î+γQ2+(1-γ)F*(s)diag⁡{fr,0(s),r=N,N-1,…,0}Q3]e.Proof.
The sojourn time of a customer who requiresm servers at Station 2 consists of:(i)
the service time of the customer when at leastm servers are available at the time of the request;(ii)
zero time, if the required number of servers is not available and the customer leaves the system;(iii)
the blocking time and the service time of a customer when the required number of servers is not available and the customer awaits the release of sufficiently many servers.
Note that we assume that the service of typem is performed by m servers independently of each other and finishes when all m servers complete the service. The distribution of this service time is defined by the Laplace-Stieltjes transform fm,0(s),m=1,N¯.
Taking into account this fact together with (i)–(iii) and using the matrix notation, we obtain expression (4.16) for the vector v2(s). In the expression, the first summand corresponds to the case (i), the second and the third summands give the Laplace-Stieltjes transform under study in the cases (ii) and (iii), respectively.


### 4.1.3. The Virtual Sojourn Time in the System

Theorem 4.5.
The Laplace-Stieltjes transform of the virtual sojourn time distribution in the system is given by(4.17)v(s)=v1(s)(IN+1⊗eW¯)v2(s),
where the vectors v1(s) and v2(s), respectively, are as given in (4.10) and (4.16).Proof.
Formula (4.17) readily follows from the structure of the virtual sojourn time in the system which consists of the virtual sojourn time at Station 1 and the sojourn time at Station 2.


## 4.1.1. The Virtual Sojourn Time at Station 1

The virtual sojourn time at Station 1 consists of (i) theresidual time from an arbitrary time instant (associated with virtual customer arrival) to the next service completion epoch at Station 1 (ii) the generalized service times of customers staying in the queue at an arbitrary time, and (iii) the generalized service time of the virtual customer.First, we study theresidual time. To this end, we consider the process χt={it,mt,rt,νt,ṽt},t≥0, whose components are defined as follows: it is the number of customers in Station 1 (including the blocked customer, if any), mt takes values 0, 1, 2, respectively, based on the server at Station 1 is idle, busy, or blocked at time t, νt is the state of the BMAP, rt is the number of busy servers at Station 2 just before the service completion epoch following the time t, ṽt is the residual time from t to that service completion epoch.Using the definition of semiregenerative processes given in [17], it can be verified that the process χt is a semi-regenerative one with the embedded Markov renewal process {ξn,tn},n≥1. Let(4.1)Ṽ(i,m,r,ν,x)=lim⁡t→∞P{it=i,mt=m,rt=r,νt=ν,ṽt<x},i≥0,m=0,2¯,r=0,N¯,ν=0,W¯,x≥0,
be the stationary distribution of the process χt,t≥0.From [17], the limits in (4.1) exist if the process {ξn,tn},n≥1, is irreducible aperiodic recurrent and the value τ of the mean inter-departure time at Station 1 (given by (3.16)) is finite. All these conditions hold if inequality (3.7) is satisfied.LetṼ(i,m,x) be the row vector of the steady state probabilities Ṽ(i,m,r,ν,x) arranged according to the lexicographic order of the components (r,ν), and let ṽ(i,m,s) be the corresponding vector of the Laplace-Stieltjes transforms, that is, let ṽ(i,m,s)=∫0∞e-sxdṼ(i,m,x),i≥0,m=0,2¯.Lemma 4.1.
The vector Laplace-Stieltjes transformsṽ(i,m,s) are calculated by
(4.2)ṽ(0,1,s)=0,ṽ(i,0,s)=0,i>0,ṽ(0,2,s)=0,(4.3)ṽ(0,0,s)=-τ-1π0Q̂(Δ⊕D0)-1[B*(s)⊗IW¯],(4.4)ṽ(i,1,s)=τ-1{π0∑k=1i[-Q̂(Δ⊕D0)-1D̃k+(1-γ)FkQ̃3]×∫0∞(eΔu⊗IW¯)∫0uIN+1⊗P(i-k,y)e-s(u-y)dydB(u)+∑j=1iπj[(Q̃1+γQ̃2)∫0∞(eΔu⊗IW¯)∫0uIN+1⊗P(i-j,y)e-s(u-y)dydB(u)+(1-γ)∑k=0i-jFkQ̃3∫0∞eΔu⊗P(i-k-j,y)e-s(u-y)dydB(u)]},(4.5)ṽ(i,2,s)=τ-1(1-γ)∑j=0i-1πj(∫0∞dF(u)⊗IW¯)Q̃3×∫0ue-s(u-y)[IN+1⊗P(i-j-1,y)]dy[B*(s)⊗IW¯],i>0.Proof.
Letκj(i,m,x,t)(r,ν;r′,ν′) denote the conditional probability that, given time 0 is an instant of the service completion at Station 1 and the embedded Markov chain ξn is in the state (j,r,ν) at that time, the next service completion epoch at Station 1 occurs later than t, the discrete components of the process χt take values (i,m,r′,ν′) at time t and the continuous-time component ṽt<x.
Let us arrange the probabilitiesκj(i,m,x,t)(r,ν;r′,ν′), for fixed values i,j,m, according to the lexicographic order of the states (r,ν;r′,ν′) and form the square matrices
(4.6)K̃j(i,m,x,t)=(κj(i,m,x,t)(r,ν;r′,ν′))ν,ν′=0,W¯;r,r′=0,N¯.
Then, using the ergodic theorem for semi-regenerative processes, (see [17, Theorem  6.12]), the probability vectors Ṽ(i,m,x) can be related to the stationary distribution πj,j≥0, of the embedded Markov chain ξn,n≥1, by
(4.7)Ṽ(i,m,x)=τ-1∑j=0∞πj∫0∞K̃j(i,m,x,t)dt,i≥0,m=0,2¯.
The corresponding vector Laplace-Stieltjes transforms ṽ(i,m,s) are defined by
(4.8)ṽ(i,m,s)=τ-1∑j=0∞πj∫0∞K̃j*(i,m,s,t)dt,i≥0,m=0,2¯,
where K̃j*(i,m,s,t)=∫0∞e-sxdK̃j(i,m,x,t).
From (4.8), formulas (4.2) follow immediately when we note that K̃j*(i,m,s,t)=0 for the range of arguments {j≥0,i=0,m=1,2} and {j≥0,i>0,m=0}.
Letm=1. The lengthy but straightforward expressions for the matrices K̃j*(i,1,s,t), i>0, are presented in Appendix A. Substituting these expressions into (4.8) and after routine algebraic manipulations including rearranging the order of integration, we get formula (4.4) for the vectors ṽ(i,1,s),i>0. Similar calculations yield (4.3) and (4.5).Further, we study thegeneralized service time distribution at Station 1.LetB̂(x) be the matrix distribution function of generalized service time. More specifically, let B̂(x)=(B̂(x)r,r′)r,r′=0,N¯, where B̂(x)r,r′=P{tn+1-tn<x,rn+1=r′∣rn=r,in≠0}. Denote ℬ(s)=∫0∞e-stdB̂(t).Lemma 4.2.
The matrix Laplace-Stieltjes transform of the generalized service time distribution at Station 1 is calculated as(4.9)B(s)=[Q1+γQ2+(1-γ)F*(s)Q3]B*(s),
where F*(s)=∫0∞e-stdF(t).Proof.
To prove, we need to analyze the structure of the generalized service time. The generalized service time of a tagged customer is just the service time of the customer at Station 1 if the previous customer did not block the server of this station. In this case, the matrix Laplace-Stieltjes transform of the generalized service time distribution is calculated by(Q1+γQ2)B*(s). However, when blocking occurs, the generalized service time consists of the time during which the server is blocked by the previous customer and the service time of the tagged customer. The corresponding Laplace-Stieltjes transform is defined by (1-γ)∫0∞e-stdF(t)Q3B*(s). The stated result (4.9) follows immediately.Now we are ready to derive the equation for the vector Laplace-Stieltjes transformv1(s) of the distribution of the virtual sojourn time at Station 1. Let v1(r,ν,x) be the probability that, at an arbitrary epoch, the BMAP is in state ν, the virtual sojourn time at Station 1 is less than x, and the number of busy servers at Station 2 just before the end of the virtual sojourn time is r. Then v1(s) is defined as a vector of Laplace-Stieltjes transforms v1(r,ν,s)=∫0∞e-sxdv1(r,ν,x) written in lexicographic order.Theorem 4.3.
The vector Laplace-Stieltjes transformv1(s) satisfies the equation
(4.10)v1(s)A(s)=π0Φ(s),
where
(4.11)A(s)=sI+∑r=0∞Br(s)⊗Dr,Φ(s)=τ-1Q̂(Δ⊕D0)-1(Δ⊗IW¯-sI)[B*(s)⊗IW¯].Proof.
As mentioned above, the virtual sojourn time at Station 1 consists of the residual time from an arbitrary timet to the next service completion epoch, the generalized service times of customers that await for a service at time t, and the generalized service time of the virtual customer.
Taking into account the structure of the virtual sojourn time and using the law of total probability, we express the vector Laplace-Stieltjes transformv1(s) as follows:(4.12)v1(s)=ṽ(0,0,s)+∑i=1∞ṽ(i,1,s)[Bi(s)⊗IW¯]+∑i=1∞ṽ(i,2,s)[Bi-1(s)⊗IW¯].
Further, we multiply (4.12) by the matrix sI+∑r=0∞ℬr(s)⊗Dr, and, after some laborious calculations, get(4.13)v1(s)(sI+∑r=0∞Br(s)⊗Dr)=τ-1{π0∑i=0∞[Ci+∑j=1i+1πjYi-j+1][Bi+1(s)⊗IW¯]+π0Q̂[-(Δ⊕D0)-1(sI+D̃0)+I][B*(s)⊗IW¯]-∑j=0∞πj[Bj+1(s)⊗IW¯]}.
Multiplying the balance equations for stationary probability vectors πi of the form
(4.14)πi=π0Ci+∑l=1i+1πlYi-l+1,i≥0,
by ℬi+1(s)⊗IW¯ and summing over i we obtain
(4.15)∑i=0∞πi[Bi+1(s)⊗IW¯]=π0∑i=0∞Ci[Bi+1(s)⊗IW¯]+∑i=0∞∑j=1i+1πjYi-j+1[Bi+1(s)⊗IW¯].
Using (4.15) to simplify equation (4.13), we obtain (4.10).


## 4.1.2. The Sojourn Time at Station 2

Letv2(s) be the column vector of the Laplace-Stieltjes transforms of the conditional sojourn time distributions at Station 2. The rth entry of this vector is the Laplace-Stieltjes transform of the sojourn time distribution of a customer at Station 2 given that the number of busy servers is equal to r just before the end of the sojourn time of this customer at Station 1.Lemma 4.4.
The vector Laplace-Stieltjes transform of the sojourn time distribution at Station 2 is given by(4.16)v2(s)=[Q4(F*(s)+I)Î+γQ2+(1-γ)F*(s)diag⁡{fr,0(s),r=N,N-1,…,0}Q3]e.Proof.
The sojourn time of a customer who requiresm servers at Station 2 consists of:(i)
the service time of the customer when at leastm servers are available at the time of the request;(ii)
zero time, if the required number of servers is not available and the customer leaves the system;(iii)
the blocking time and the service time of a customer when the required number of servers is not available and the customer awaits the release of sufficiently many servers.
Note that we assume that the service of typem is performed by m servers independently of each other and finishes when all m servers complete the service. The distribution of this service time is defined by the Laplace-Stieltjes transform fm,0(s),m=1,N¯.
Taking into account this fact together with (i)–(iii) and using the matrix notation, we obtain expression (4.16) for the vector v2(s). In the expression, the first summand corresponds to the case (i), the second and the third summands give the Laplace-Stieltjes transform under study in the cases (ii) and (iii), respectively.


## 4.1.3. The Virtual Sojourn Time in the System

Theorem 4.5.
The Laplace-Stieltjes transform of the virtual sojourn time distribution in the system is given by(4.17)v(s)=v1(s)(IN+1⊗eW¯)v2(s),
where the vectors v1(s) and v2(s), respectively, are as given in (4.10) and (4.16).Proof.
Formula (4.17) readily follows from the structure of the virtual sojourn time in the system which consists of the virtual sojourn time at Station 1 and the sojourn time at Station 2.


## 4.2. The Actual Sojourn Time

Letv1(a)(s) and v(a)(s) be the Laplace-Stieltjes transforms of the distribution of the actual sojourn time at Station 1 and in the whole system.Theorem 4.6.
The Laplace-Stieltjes transform of the actual sojourn time distribution at Station 1 is calculated as follows:(4.18)v1(a)(s)=λ-1v1(s)∑k=0∞[Bk(s)(B(s)-I)-1⊗Dk]e.Proof.
The actual sojourn time at Station 1 of an arbitrary-tagged customer, who arrived in a group of sizek and placed at the jth position within the group, consists of (a) the actual sojourn time at Station 1 of the first customer in the group, which coincides with the virtual sojourn time at Station 1; (b) the generalized service times at Station 1 of the j-2 customers of the group who arrived with the tagged customer; (c) the generalized service time of the tagged customer at Station 1.
The vector Laplace-Stieltjes transform of the sojourn time distribution at Station 1 of the first customer of thek-size group that contains the tagged customer is evidently given by the vector v1(s)(IN+1⊗kDke/λ).
Assuming that an arbitrary customer arriving in a group of sizek is placed on the jth position with probability 1/k and using the law of total probability, we immediately obtain the following expression:
(4.19)v1(a)(s)=∑k=1∞v1(s)(IN+1⊗kDkeλ)∑j=1k1kBj-1(s)e.
After some algebraic manipulations (4.19) is reduced to (4.18).Corollary 4.7.
The Laplace-Stieltjes transform of the actual sojourn time distribution in the whole system is calculated as follows:(4.20)v(a)(s)=λ-1v1(s)∑k=0∞[Bk(s)(B(s)-I)-1⊗Dk](IN+1⊗eW¯)v2(s).


## 4.3. Moments of the Sojourn Time Distribution

The formulas for the moments of the virtual sojourn time distribution can be obtained by differentiating the expression in (4.10) at the point s=0. This requires the calculation of the derivatives of v1(s) at the point s=0. However, the matrix A(s) in (4.10) is singular at the point s=0 and calculation of the derivatives at this point is the nontrivial task. The results given below allow one to develop a procedure for calculating the required derivatives.We will use notationv1(m)(s) for the mth derivative of the vector v1(s),m≥1,and set v1(0)(s)=v1(s). Similar notations will be used for other functions of s.Theorem 4.8.
Let∫0∞tmdB(t)<∞,m=1,M+1¯, where M is an arbitrary positive integer. Then the vectors v1(m)(0),m=1,M¯, are computed recursively by
(4.21)v1(m)(0)=[(π0Φ(m)(0)-∑l=0m-1(ml)v1(l)(0)A(m-l)(0))Ĩ+1m+1(π0Φ(m+1)(0)-∑l=0m-1(m+1l)v1(l)(0)A(m+1-l)(0))eê]Ã-1,
with the initial condition
(4.22)v1(0)(0)=v1(0)=[τ-1π0Q̂(Δ⊕D0)-1(ΔB*(0)⊗IW¯)Ĩ+pidleê]Ã-1,
where Ã=A(0)Ĩ+A′(0)eê.Proof of the theorem is presented in AppendixB. In what follows we assume that ∫0∞tkdB(t)<∞,k=1,2¯.Corollary 4.9.
The mean virtual sojourn time at Station 1 is given by(4.23)v¯1={[τ-1π0Q̂(Δ⊕D0)-1((B*(0)-ΔB*′(0))⊗IW¯)+v1(0)A′(0)]Ĩ+[pidleb1+12v1(0)A′′(0)e]ê}Ã-1e,
where the vector v1(0) is defined by formula (4.22).Proof follows from the relationv¯1=-v1′(0)e and formula (4.21).Theorem 4.10.
The mean virtual sojourn time in the system is given by(4.24)v¯=v¯1+v1(0)(IN+1⊗eW¯)v¯2,
where v1(0) and v¯1 are given in (4.22) and (4.23), respectively, and v¯2 is a vector of conditional means of the sojourn time at Station 2,
(4.25)v¯2=-{Q4F(1)êT+(1-γ)[F(1)-Ediag⁡{∑l=1r(lμ)-1,r=N,N-1,…,0}]Q3e},
where F(1) is given by formula (C.3) below.Proof.
To calculate the valuev¯ we differentiate the expression in (4.17). Setting s=0 and replacing the sign, we have that
(4.26)v¯=-v1′(0)(IN+1⊗eW¯)v2(0)-v1(0)(IN+1⊗eW¯)v2′(0).
Putting s=0 in (4.16) we get v2(0)=[Q4(E+I)Î+γQ2+(1-γ)EQ3]e=e. This implies that the first term in the right-hand side of (4.26) is equal to -v1′(0)e=v¯1. Using the relation v¯2=-v2′(0) and differentiating (4.16) at the point s=0, we readily verify that v¯2 has form (4.25). This completes the proof.Theorem 4.11.
The mean actual sojourn time at Station 1 is given by(4.27)v¯1(a)=-λ-1{v1′(0)(e⊗∑k=1∞kDke)+v1(0)∑k=1∞(IN+1⊗Dke)[∑n=1k-1∑l=0n-1Bl(0)B′(0)e]}.Proof of the theorem follows from the relationv¯1(a)=-dv1(a)(s)/ds|s=0 and formula (4.19).Corollary 4.12.
The mean actual sojourn time in the system is given by(4.28)v¯a=v¯1(a)+λ-1v1(0)∑k=1∞(IN+1⊗Dke)∑n=0k-1Bn(0)v¯2.


## 5. Numerical Examples

In this section, we demonstrate feasibility of the algorithms developed here and show numerically some interesting features of the system under consideration.Experiment 1.
In this experiment, we investigate the impact of coefficient of variation in the service process at Station 1 on the main performance measures of the system.To this end, we consider four service processes with the same mean service timeb1=0.1, but different values for the coefficient of variation, cvar. The first process is coded as D and corresponds to the deterministic service time distribution. The second process is coded as M and corresponds to the exponential service time. The third and the fourth service processes are coded as HM2(1), HM2(2) and correspond to hyperexponential service time distributions of order 2. These distributions are defined by the mixing probabilities (0.05,0.95) and the intensities 0.62025,48.9998 in the case of HM2(1) and (0.98,0.02) and the intensities 10000,0.2 in the case of HM2(2). The coefficients of variation of processes D, M, HM2(1), and HM2(2) are, respectively, equal to 0,1,5,9.95.The input process is defined by the matrices(5.1)D0=(-15.73270.60620.59240.5178-2.28970.46790.59710.5653-1.9597),D=(14.15020.30210.08180.10711.0320.16460.08580.19790.5136).
The matrices Dk,k=1,5¯, are calculated as follows: Dk=Dhk-1(1-h)/(1-h5), where h=0.8. Then we normalize the matrices Dk,k=0,5¯, so as to get the arrival rate λ=1. This BMAP has coefficient of correlation ccor=0.2. The other parameters of the system are as follows: N=5,μ=0.8,γ=0.5,q0=0.1,q1=q2=0.3,q3=q4=q5=0.1.We vary the mean service timeb1 for all considered service processes in the interval [0.1,0.95] by scaling appropriately. The coefficients of variation do not change under such scaling. Note that, as it was mentioned in Remark 3.5 above, the system load ρ depends not only on the mean service time, but also on the variance of the service time. In Table 1, the value of ρ is given as function of b1 and cvar. The values of ρ that exceed 1 are printed in bold face. The tandem queueing system is not stable for these values.Figures1 and 2 show the dependence of the main performance measures of the system on the value of the mean service time b1 for service processes with different service time distribution. From these figures, it is very clear that the key performance measures of the system are very sensitive with respect to the service time variance. We also ran other examples, besides the one presented here, involving Erlangian and uniform distributions for the service time distribution. Since all these distributions have coefficient of variation in the range (0,1), the corresponding curves, as expected, were located between the two lower curves in Figures 1 and 2.The mean virtual and actual sojourn time as functions of the mean service time for different service time distributions.
(a)(b)The loss probability and the mean number of busy servers at Station 2 as functions of the mean service time for different service time distributions.
(a)(b)Experiment 2.
Here we solve numerically the following optimization problem. Find an optimal choice for the numberN of servers at Station 2 that will minimize the expected total cost per unit of time:
(5.2)J=J(N)=aN+c1λPloss+c2v¯(a),
where a is the cost of utilization per unit time of a server at Station 2 (maintenance cost), c1 is the cost of a customer leaving the system after a service at Station 1 due to lack of required servers, and c2 is the cost per unit of time of holding (the sojourn time) an arbitrary customer in the system (holding cost).Using theMAP characterized by the matrices(5.3)D0=(-6.745385.45412×10-65.45412×10-6-0.219455),D=(6.7006850.0446950.1224270.097023),
we construct the BMAP with the matrices Dk,k=0,5¯, similar to Experiment 1 and normalize so as to have λ=3. The BMAP has the coefficient of correlation ccor=0.2 and the coefficient of variation cvar=3.5.The service time distribution at Station 1 is assumed to be Erlang of order 3 with parameter 20. The probabilityγ is taken to be 0.5. The components of the vector q are: q0=0.1,q1=0.9,qm=0,m=2,N¯. The various costs are taken as follows: a=5,c1=50,c2=3.The objective function,J, as a function of the number N of servers under different service rates μ is plotted in Figure 3. Table 2 contains the values of the objective function.Table 2
The value of the objective function for different number of servers and different service rate at Station 2.
μ=1μ=2μ=3μ=4μ=5N=1∞457.034674.935151.031142.5146N=2353.737848.634934.113027.082622.9120N=371.048034.656525.844122.265620.6657N=451.782230.406825.820224.600824.2393N=545.255631.208029.392529.114129.0630N=643.075034.643734.098034.053534.0486N=743.686739.071539.053339.050439.0476N=847.256444.052744.051444.049544.0474N=951.985449.052149.051149.049149.0473N=1057.005655.045854.600154.048954.0472Figure 3
The objective function as a function of the number of servers in Station 2 for different service rates.The optimal valuesJ* of the objective function for each of the five service rates are displayed in bold face. It is seen from Figure 3 and Table 2 that as the service rate decreases from 5 to 1, the optimal number of servers N* increases from 3 to 6. The relative gain of the optimal configuration in comparison to a system with an arbitrary number N of the servers at Station 2 is defined as Rrel⁡(N)=((J(N)-J*)/J*)100%.We now focus on the result of optimal value in the caseμ=5. It is seen from Figure 3 and Table 2 that the optimal value of the objective function J* is 20.6657 and the optimal number of servers N*=3. It should be also noted that for the case under consideration the minimal relative gain is more than 10% if we install the optimal number of servers N*=3 instead of 2 servers and maximal relative gain is more than 161% if we use N*=3 servers instead of 10 servers.Experiment 3.
In this experiment, we show that the correlation in the input flow has a great impact on the performance measures of the system. In addition to the BMAP defined in the first experiment and having the coefficient of correlationccor=0.2, let us consider two another BMAPs, having the same mean arrival rate, but different coefficients of correlation. These BMAPs are defined by the matrices D0 and D1=D, from which the matrices Dk,k=0,5¯, are defined in the same way as in Experiment 1.TheBMAP having the coefficient of correlation ccor=0.1 is characterized by the matrices(5.4)D0=(-13.33460.58860.61730.6927-2.44660.42290.68230.4144-1.6354),D=(11.54690.36310.21870.38420.86590.08090.28520.04250.2111).
The BMAP having the coefficient of correlation ccor=0.3 is defined by the matrices(5.5)D0=(-25.53980.39330.36120.1452-2.23220.20000.29600.3874-1.7526),D=(24.24210.46690.07630.03411.66680.18610.00900.25550.8047).
All the three defined BMAPs have the coefficient of variation cvar=2.In addition, we consider theBMAP which is a group Poisson process. It has the same mean arrival rate as three other BMAPs, coefficient of correlation ccor=0 and the coefficient of variation cvar=1.The service time at Station 1 has the Erlangian distribution of order 3 with parameter 20. The mean service timeb1=3/20 and the squared coefficient of variation cvar2=1/3.The number of servers at Station 2N=5. Service rate μ is equal to 5. Probability γ that a customer will await the release of servers is equal to 0.5. The components of the vector q, which defines the type of the customer, are q0=0.1,q1=q2=0.3, and q3=q4=q5=0.1.Figures4 and 5 illustrate the dependence of the mean virtual sojourn time v¯, the mean actual sojourn time v¯a, the loss probability Ploss, and the mean number Nbusy of busy servers at Station 2 on the system load ρ. The load ρ varies by means of scaling the fundamental rate λ. Note that the coefficients of correlation and variation of the BMAP do not change under such scaling.The mean virtual and actual sojourn time as functions of the system load for theBMAPs with different correlation.
(a)(b)The loss probability and the mean number of busy servers at Station 2 as functions of the system load for theBMAPs with different correlation.
(a)(b)Figures4 and 5 confirm the fact that values of v¯, v¯a, Ploss, and Nbusy increase when the system load, ρ, increases. We also note that, under the same scenario for the system load the correlation of the interarrival times shows a strong (negative) impact on the system performance characteristics.


## 6. Conclusion

In this paper, theBMAP/G/1→•/M/N/0 tandem queue with heterogeneous customers is studied. The system is studied by looking at selected embedded epochs. The condition for the existence of the stationary distribution is derived. Expressions for the loss probability, blocking probability, and some other performance characteristics of the system are obtained. The Laplace-Stieltjes transforms of the distribution of the virtual and the actual sojourn time at both stations as well as at the whole system are derived. Although the required analytical derivations are very complicated and cumbersome, the resulting formulas have very simple forms. The procedure for calculating the moments of the virtual and the actual sojourn time distribution is elaborated. Illustrative numerical results highlight the important role played by the variance of the service time. An optimization problem to illustrate the usefulness of such problems in practice involving tandem queues is discussed. The results of this paper can be applied to areas such as capacity planning, performance evaluation, and optimization of real-world tandem queues and two-node networks.



---

*Source: 324604-2011-12-15.xml*