A Tandem BMAP/𝐺/1→•/𝑀/𝑁/0 Queue with Group Occupation of Servers at the Second Station

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.

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 the BMAP input were considered in [6–10]. The papers [6, 7, 10] are devoted to the MAP/PH/1→•/𝐺/1 system with blocking. In [8], the tandem queues BMAP/𝐺/1/𝑁→•/PH/1/𝑀−1 with losses are studied. The tandem-queue of the BMAP/𝐺/1→•/PH/1/𝑀−1 type with losses and feedback has been studied in [9].

In the present paper, we consider a BMAP/𝐺/1→•/𝑀/𝑁/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 with MAP 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 Section 2, 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 the BMAP/𝐺/1 queue. That is, the arrivals to Station 1 are described by a BMAP. The BMAP is defined by the underlying process 𝜈_𝑡,𝑡≥0, which is an irreducible continuous time Markov chain with state space {0,…,𝑊}, and with the matrix generating function 𝐷(𝑧)=∑^∞_𝑘=0𝐷_𝑘𝑧^𝑘,|𝑧|≤1. Arrivals occur only at epochs of the jumps in the underlying process 𝜈_𝑡,𝑡≥0. The intensities of the transitions of the process 𝜈_𝑡 accompanied by a batch of size 𝑘 are defined by the matrices 𝐷_𝑘,𝑘≥0. The matrix 𝐷(1) is the infinitesimal generator of the process 𝜈_𝑡. The stationary distribution vector 𝜽 of this process satisfies the equations 𝜽𝐷(1)=𝟎,𝜽𝐞=1, where 𝐞 is a column vector consisting of 1′s, and 𝟎 is a row vector of 0^’s.

The average intensity 𝜆 (fundamental rate) of the BMAP is given by 𝜆=𝜽𝐷^′(𝑧)|_𝑧=1𝐞. We assume that 𝜆<∞. The average intensity 𝜆_𝑏 of group arrivals is defined by 𝜆_𝑏=𝜽(−𝐷₀)𝐞. The coefficient of variation, 𝑐_var, of intervals between successive group arrivals is defined by 𝑐²_var=2𝜆_𝑏𝜽(−𝐷₀)⁻¹𝐞−1. The coefficient of correlation 𝑐_cor of the successive intervals between group arrivals is given by 𝑐_cor=(𝜆_𝑏𝜽(−𝐷₀)⁻¹(𝐷(1)−𝐷₀)(−𝐷₀)⁻¹𝐞−1)/𝑐²_var. 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 distribution 𝐵(𝑡), Laplace-Stieltjes transform 𝛽(𝑠)=∫^∞₀𝑒^−𝑠𝑡𝑑𝐵(𝑡), and finite first moment 𝑏₁=∫^∞₀𝑡𝑑𝐵(𝑡).

After receiving a service at Station 1, the customer proceeds to Station 2. At this station, there are 𝑁 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 𝑞_𝑚,𝑞_𝑚≥0,𝑚=0,𝑁,∑^𝑁_𝑚=0𝑞_𝑚=1, the customer will require exactly 𝑚 servers to provide a service at Station 2 and will be called type 𝑚 customer. Here and in the sequel, notation such as 𝑚=0,𝑁, means that 𝑚 assumes values from the set {0,1,…,𝑁}. 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 𝑞₀≠1. Otherwise, the queue under consideration will be reduced to the BMAP/𝐺/1 queue which has been studied extensively.

If the customer is of type 𝑚, 𝑚=1,𝑁, and the required number of servers is available, the customer’s service will begin immediately. Each of these 𝑚 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)𝐼 is an identity matrix of appropriate dimension;(ii)⊗ and ⊕ are symbols of the Kronecker product and sum of matrices;(iii)𝐷_𝑘=𝐼_𝑁+1⊗𝐷_𝑘,𝑘≥0,𝐷(𝑧)=∑^∞_𝑘=0𝐷_𝑘𝑧^𝑘,|𝑧|≤1; (iv)𝑃(𝑗,𝑡),𝑗≥0, is a matrix function defined by the expansion ∑^∞_𝑗=0𝑃(𝑗,𝑡)𝑧^𝑗=𝑒^{𝐷(𝑧)𝑡};(v)𝐹(𝑡)=(𝐹_𝑟,𝑟′(𝑡))_{𝑟,𝑟′=0,𝑁}, where 𝐹_𝑟,𝑟′(𝑡)=0 for 𝑟≤𝑟^′ and for 𝑟>𝑟^′,𝐹_𝑟,𝑟′(𝑡) is the generalized Erlang distribution function with the Laplace-Stieltjes transform 𝑓_𝑟,𝑟′(𝑠)=∏^𝑟_{𝑙=𝑟′+1}𝑙𝜇(𝑙𝜇+𝑠)⁻¹;(vi)𝑄_𝑚,𝑚=1,4, are square matrices: 𝑄₂=diag{∑^𝑁_{𝑚=𝑁−𝑟+1}𝑞_𝑚,𝑟=0,𝑁}, 𝑄₁=𝑞₀𝑞₁…𝑞_𝑁0𝑞₀…𝑞_𝑁−1⋮⋮⋱⋮00…𝑞₀,𝑄₃=0…0𝑞_𝑁0…0𝑞_𝑁−1⋮⋱⋮⋮0…0𝑞₀,𝑄₄=𝑞₀𝑞₁…𝑞_𝑁−1𝑞_𝑁𝑞₀𝑞₁…𝑞_𝑁−10𝑞⋮⋮⋱⋮⋮₀;0…00 (2.1)(vii)𝑄_𝑚=𝑄_𝑚⊗𝐼_𝑊,𝑚=1,3,𝑊=𝑊+1; (viii)𝑄=𝑄₁+𝛾𝑄₂+(1−𝛾)∫^∞₀(𝑑𝐹(𝑡)⊗𝑒^𝐷0𝑡)𝑄₃,𝑄=𝑄₁+𝛾𝑄₂+(1−𝛾)𝐸𝑄₃; (ix)𝐸=,̂,̃̂̂().00⋯0010…00⋮⋮⋱⋮⋮11⋯10𝐼=10⋯000⋯0⋮⋮⋱⋮00⋯0𝐼=𝐼−𝐼,𝐞=1,0,…,0 (2.2)

3. The Stationary Distribution of the Embedded Markov Chain

Let 𝑡_𝑛 denote the time of the 𝑛th service completion at Station 1. Consider the process 𝜉_𝑛={𝑖_𝑛,𝑟_𝑛,𝜈_𝑛}, 𝑛≥1, where 𝑖_𝑛,𝑖_𝑛≥0, is the number of customers at Station 1 (not counting the blocked customer, if any) at epoch 𝑡_𝑛+0; 𝑟_𝑛,𝑟_𝑛=0,𝑁, is the number of busy servers at Station 2 at epoch 𝑡_𝑛−0; 𝜈_𝑛,𝜈_𝑛=0,𝑊, is the state of the BMAP at epoch 𝑡_𝑛.

It is easy to verify that the process 𝜉_𝑛={𝑖_𝑛,𝑟_𝑛,𝜈_𝑛}, 𝑛≥1, is a Markov chain. Enumerating the states of this Markov chain in lexicographic order, and denoting by 𝑃_𝑙,𝑘,𝑙,𝑘≥0, the square matrix of order (𝑊+1)(𝑁+1) governing the transition probabilities of the chain from the set of states {𝑙,⋅,⋅} to the set {𝑘,⋅,⋅}, the following lemma gives the entries of the transition probability matrix of the Markov chain 𝜉_𝑛.

Lemma 3.1. The transition probability matrix of the chain 𝜉_𝑛,𝑛≥1, has the following block structure: 𝑃=𝑃_{𝑙,𝑘𝑙,𝑘≥0}=𝐶₀𝐶₁𝐶₂𝐶₃⋯Y₀𝑌₁𝑌₂𝑌₃⋯0𝑌₀𝑌₁𝑌₂⋯00𝑌₀𝑌₁⋯,⋮⋮⋮⋮⋱ (3.1) where 𝐶_𝑖=^𝑖+1_𝑘=1−𝑄Δ⊕𝐷₀⁻¹𝐷_𝑘+1−𝛾𝐹_𝑘𝑄₃Ω_{𝑖−𝑘+1},𝑌_𝑖=𝑄₁+𝛾𝑄₂Ω_𝑖+1−𝛾^𝑖_𝑘=0𝐹_𝑘𝑄₃Ω_𝑖−𝑘,Ω_𝑗=^∞₀𝑒^Δ𝑡⊗𝑃𝑗,𝑡𝑑𝐵(𝑡),𝐹_𝑗=^∞₀(𝑡).𝑑𝐹⊗𝑃𝑗,𝑡,𝑗≥0,Δ=000⋯00𝜇−𝜇0⋯0002𝜇−2𝜇⋯00⋮⋮⋮⋱⋮⋮000⋯𝑁𝜇−𝑁𝜇 (3.2)

Proof. First, we write the transition probability matrices 𝐶_𝑖,𝑌_𝑖 in block forms as 𝐶_𝑖=(𝐶^{(𝑟,𝑟′)}_𝑖)_{𝑟,𝑟′=0,𝑁},𝑌_𝑖=(𝑌^{(𝑟,𝑟′)}_𝑖)_{𝑟,𝑟′=0,𝑁}, where the blocks 𝐶^{(𝑟,𝑟′)}_𝑖, 𝑌^{(𝑟,𝑟′)}_𝑖 correspond to the transitions of the number of busy servers from 𝑟 to 𝑟^′ at Station 2.
Denote by 𝛿_𝑟,𝑟′(𝑡) the probability that during the time interval of the length 𝑡 the number of busy servers at Station 2 decreases from 𝑟 to 𝑟^′ 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 𝑃(𝑗,𝑡) gives the probability that 𝑗 customers arrive in the BMAP during the interval (0,𝑡] and the state of the BMAP at epoch 𝑡 is 𝜈′ given 𝜈₀=𝜈.
The (𝜈,𝜈^′)th entry of the matrix ∫^∞₀𝛿_𝑟,𝑟′(𝑡)𝑃(𝑗,𝑡)𝑑𝐵(𝑡) gives the probability that during the service time of a customer at Station 1, exactly 𝑗 customers arrive, the number of busy servers at Station 2 decreases from 𝑟 to 𝑟^′ and the BMAP has changed from 𝜈 to 𝜈^′.
The (𝜈,𝜈^′)th entry of the matrix ∫^∞₀𝛿_𝑟,𝑟′(𝑡)𝑒^𝐷0𝑡𝐷_𝑘𝑑𝑡 gives the probability that at an arbitrary time instant with the number of busy servers at Station 2 equal to 𝑟, and the BMAP in state 𝜈, the first batch to arrive is of size 𝑘; soon after that instant, the number of busy servers at Station 2 equals 𝑟′ and the BMAP is in state 𝜈′.
𝐹_𝑟,𝑟′(𝑡) is the distribution function of the time interval during which the number of busy servers at Station 2 decreases from 𝑟 to 𝑟^′ conditioned on the fact that none arrive to this station. Then the (𝜈,𝜈^′)th entry of the matrix ∫^∞₀𝑃(𝑗,𝑡)𝑑𝐹_𝑟,𝑟′(𝑡) defines the probability that exactly 𝑗 customers arrive with the BMAP moving from 𝜈 to 𝜈^′ and that the number of busy servers at Station 2 decreases from 𝑟 to 𝑟^′ during that time interval.
From the above probabilistic interpretations, analyzing the one-step transitions of the chain 𝜉_𝑛 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 𝐶^{(𝑟,𝑟′)}_𝑖,𝑌^{(𝑟,𝑟′)}_𝑖,𝑖≥0: 𝐶^{(𝑟,𝑟′)}_𝑖=^𝑁−𝑟_𝑚=0𝑞_𝑚^𝑟+𝑚_𝑙=𝑟′^∞₀𝛿_{𝑟+𝑚,𝑙}(𝑡)𝑒^𝐷0𝑡𝑑𝑡^𝑖+1_𝑘=1𝐷_𝑘^∞₀𝛿_𝑙,𝑟′(𝑡)𝑃(𝑖−𝑘+1,𝑡)𝑑𝐵(𝑡)+𝛾^𝑁_{𝑚=𝑁−𝑟+1}𝑞_𝑚^𝑟_𝑙=𝑟′^∞₀𝛿_𝑟,𝑙(𝑡)𝑒^𝐷0𝑡𝑑𝑡^𝑖+1_𝑘=1𝐷_𝑘^∞₀𝛿_𝑙,𝑟′(𝑡)𝑃()(𝑡)+𝑖−𝑘+1,𝑡𝑑𝐵1−𝛾^𝑁_{𝑚=𝑁−𝑟+1}𝑞_𝑚^∞₀𝑒^𝐷0𝑡𝑑𝐹_{𝑟,𝑁−𝑚}(𝑡)^𝑁_𝑙=𝑟′^∞₀𝛿_𝑁,𝑙(𝑡)𝑒^𝐷0𝑡×𝑑𝑡_𝑖+1𝑘=1𝐷_𝑘^∞₀𝛿_𝑙,𝑟′(𝑡)𝑃(𝑖−𝑘+1,𝑡)𝑑𝐵(𝑡)+^𝑖+1_𝑘=1^∞₀𝑃(𝑘,𝑡)𝑑𝐹_{𝑟,𝑁−𝑚}(𝑡)×^∞₀𝛿_𝑁,𝑟′(𝑡)𝑃()(𝑡,𝑌𝑖−𝑘+1,𝑡𝑑𝐵)_{(𝑟,𝑟′)𝑖}=^𝑁−𝑟_𝑚=0𝑞_𝑚^∞₀𝑃(𝑖,𝑡)𝛿_{𝑟+𝑚,𝑟′}(𝑡)(𝑡)+𝑑𝐵_{𝑁𝑚=𝑁−𝑟+1}𝑞_𝑚𝛾^∞₀𝛿_𝑟,𝑟′(𝑡)𝑃()(𝑡)+𝑖,𝑡𝑑𝐵1−𝛾^𝑖_𝑘=0^∞₀𝑃(𝑘,𝑡)𝑑𝐹_{𝑟,𝑁−𝑚}(𝑡)^∞₀𝛿_𝑁,𝑟′(𝑡)𝑃(𝑖−𝑘,𝑡)𝑑𝐵(𝑡). (3.3) In order to arrive at equations (3.2) from (3.3), we use the matrix notations introduced above and the relation: (𝛿_𝑟,𝑟′(𝑡))_{𝑟,𝑟′=0,𝑁}=𝑒^Δ𝑡, which follows from the fact that under the case when none arrive to Station 2, the process 𝑟_𝑡 governing the number of busy servers at this station is Markovian with generator Δ.

It is easy to see that the Markov chain 𝜉_𝑛 belongs to the class of 𝑀/𝐺/1 type Markov chains, see [16]. We can use this fact to derive the ergodicity condition and calculate the stationary distribution of the chain.

Let 𝐶(𝑧)=∑^∞_𝑖=0𝐶_𝑖𝑧^𝑖,𝑌(𝑧)=∑^∞_𝑖=0𝑌_𝑖𝑧^𝑖,|𝑧|≤1, be the generating functions of the transition probability matrices 𝐶_𝑖 and 𝑌_𝑖,𝑖≥0.

Corollary 3.2. The matrix generating functions 𝐶(𝑧),𝑌(𝑧) can be written as 𝐶(𝑧)=1𝑧−𝑄Δ⊕𝐷₀⁻¹𝐷(𝑧)−𝐷₀+1−𝛾𝐹(𝑧)−𝐹₀𝑄₃Ω(𝑧),𝑌(𝑧)=𝑄₁+𝛾𝑄₂+1−𝛾𝐹(𝑧)𝑄₃Ω(𝑧), (3.4)(3.5) where Ω(𝑧)=^∞_𝑛=0Ω_𝑛𝑧^𝑛=^∞₀𝑒^Δ𝑡⊗𝑒^{𝐷(𝑧)𝑡}𝑑𝐵(𝑡),𝐹(𝑧)=^∞_𝑛=0𝐹_𝑛𝑧^𝑛=^∞₀𝑑𝐹(𝑡)⊗𝑒^{𝐷(𝑧)𝑡}. (3.6)

Theorem 3.3. The necessary and sufficient condition for ergodicity of the Markov chain 𝜉_𝑛,𝑛≥1, is the fulfillment of the inequality 𝜌=𝜆𝑏₁+1−𝛾^𝑁_𝑟=1𝜗_𝑟^𝑁_{𝑚=𝑁−𝑟+1}𝑞_𝑚^𝑟_{𝑙=𝑁−𝑚+1}𝑙𝜇⁻¹<1. (3.7) Here (𝜗₁,…,𝜗_𝑁) is a part of the vector 𝝑=(𝜗₀,…,𝜗_𝑁), which is the unique solution to the system 𝝑𝑄𝐵^∗(0)=𝝑,𝝑𝐞=1, (3.8) where 𝐵^∗(𝑠)=∫^∞₀𝑒^−𝑠𝑡𝑒^Δ𝑡𝑑𝐵(𝑡).

Proof. It can be verified that the matrix 𝑌(1) is irreducible. Hence, from [16], the necessary and sufficient condition for ergodicity of the chain 𝜉_𝑛 is the fulfillment of the inequality 𝐱𝑌^′(1)𝐞<1, (3.9) where the vector 𝐱 is the unique solution of the system 𝐱𝑌(1)=𝐱,𝐱𝐞=1. (3.10)
The theorem will be proven if we show that inequality (3.9) is equivalent to inequality (3.7).
Let the vector 𝐱 be of the form 𝐱=𝝑⊗𝜽. (3.11) By the direct substitution into the system (3.10), where 𝑌(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 𝑧=1 and substituting the resulting expression for 𝑌^′(1) and the vector 𝐱 of form (3.11) into the inequality (3.9), we get 𝜌=𝜆𝑏₁+1−𝛾𝝑^∞₀𝑡𝑑𝐹(𝑡)𝑄₃𝐞<1. (3.12) The stated expression in (3.7) follows from (3.12) and the expression for ∫^∞₀𝑡𝑑𝐹(𝑡) 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−𝛾)∑^𝑁_𝑟=1𝜗_𝑟∑^𝑁_{𝑚=𝑁−𝑟+1}𝑞_𝑚∑^𝑟_{𝑙=𝑁−𝑚+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 𝑏₁ 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.

In what follows, we assume that the inequality (3.7) holds true.

Denote the stationary state probabilities of the Markov chain 𝜉_𝑛={𝑖_𝑛,𝑟_𝑛,𝜈_𝑛} by 𝜋(𝑖,𝑟,𝜈),𝑖≥0,𝑟=0,𝑁,𝜈=0,𝑊. Introduce the notation for the row vectors of these probabilities𝝅(𝑖,𝑟)=(𝜋(𝑖,𝑟,0),𝜋(𝑖,𝑟,1),…,𝜋(𝑖,𝑟,𝑊)),𝝅_𝑖=(𝝅(𝑖,0),𝝅(𝑖,1),…,𝝅(𝑖,𝑁)),𝑖≥0. (3.13) Let also 𝚷(𝑧)=∑^∞_𝑖=0𝝅_𝑖𝑧^𝑖,|𝑧|≤1, be the vector generating function of vectors 𝝅_𝑖,𝑖≥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 epochs 𝐿=𝚷^′(1)𝐞.(ii)The vector of the stationary distribution of the number of busy servers at Station 2 at the service completion epoch at Station 1: 𝐫=𝚷(1)(𝐼_𝑁+1⊗𝐞_𝑊).(iii)The mean number of busy servers at Station 2 at the service completion epoch at Station 1 𝑁_busy=𝐫diag𝑟,𝑟=0,𝑁𝐞. (3.14)(iv) The probability that an arbitrary customer leaves the system or causes the blocking of the server at Station 1 𝑃_loss=𝛾𝚷(1)𝑄₂𝐞,𝑃_block=1−𝛾𝚷(1)𝑄₂𝐞. (3.15)(v)The probability that the server of Station 1 is idle at an arbitrary time 𝑝_idle=𝜏⁻¹𝝅₀𝑄(−𝐷₀)⁻¹𝐞, where 𝜏 is the mean interdeparture time at Station 1, 𝜏=𝑏₁+𝝅₀𝑄−𝐷₀⁻¹𝐞−1−𝛾𝚷(1)𝐼_𝑁+1⊗𝐞𝐹⁽¹⁾𝑄₃𝐞, (3.16) the matrix 𝐹⁽¹⁾ is defined by formula (C.3) below.(vi)The probability that the server of Station 1 processes a customer at an arbitrary time 𝑝_serve=𝜏⁻¹𝑏₁.(vii)The probability that the server of Station 1 is blocked at an arbitrary time 𝑝_block=1−𝑝_idle−𝑝_serve.

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 to FIFO (first-in-first-out) discipline.

For use in the sequel, we define the generalized 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) the residual 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 the residual time. To this end, we consider the process 𝜒_𝑡={𝑖_𝑡,𝑚_𝑡,𝑟_𝑡,𝜈_𝑡,𝑣_𝑡},𝑡≥0, whose components are defined as follows: 𝑖_𝑡 is the number of customers in Station 1 (including the blocked customer, if any), 𝑚_𝑡 takes values 0, 1, 2, respectively, based on the server at Station 1 is idle, busy, or blocked at time 𝑡, 𝜈_𝑡 is the state of the BMAP, 𝑟_𝑡 is the number of busy servers at Station 2 just before the service completion epoch following the time 𝑡, 𝑣_𝑡 is the residual time from 𝑡 to that service completion epoch.

Using the definition of semiregenerative processes given in [17], it can be verified that the process 𝜒_𝑡 is a semi-regenerative one with the embedded Markov renewal process {𝜉_𝑛,𝑡_𝑛},𝑛≥1. Let𝑉(𝑖,𝑚,𝑟,𝜈,𝑥)=lim_𝑡→∞𝑃𝑖_𝑡=𝑖,𝑚_𝑡=𝑚,𝑟_𝑡=𝑟,𝜈_𝑡=𝜈,𝑣_𝑡<𝑥,𝑖≥0,𝑚=0,2,𝑟=0,𝑁,𝜈=0,𝑊,𝑥≥0, (4.1) be the stationary distribution of the process 𝜒_𝑡,𝑡≥0.

From [17], the limits in (4.1) exist if the process {𝜉_𝑛,𝑡_𝑛},𝑛≥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.

Let 𝐕(𝑖,𝑚,𝑥) be the row vector of the steady state probabilities 𝑉(𝑖,𝑚,𝑟,𝜈,𝑥) arranged according to the lexicographic order of the components (𝑟,𝜈), and let 𝐯(𝑖,𝑚,𝑠) be the corresponding vector of the Laplace-Stieltjes transforms, that is, let 𝐯(𝑖,𝑚,𝑠)=∫^∞₀𝑒^−𝑠𝑥𝑑𝐕(𝑖,𝑚,𝑥),𝑖≥0,𝑚=0,2.

Lemma 4.1. The vector Laplace-Stieltjes transforms 𝐯(𝑖,𝑚,𝑠) are calculated by 𝐯(0,1,𝑠)=𝟎,𝐯(𝑖,0,𝑠)=𝟎,𝑖>0,𝐯(0,2,𝑠)=𝟎,𝐯(0,0,𝑠)=−𝜏⁻¹𝝅₀𝑄Δ⊕𝐷₀⁻¹𝐵^∗(𝑠)⊗𝐼_𝑊,𝐯(𝑖,1,𝑠)=𝜏⁻¹𝝅₀^𝑖_𝑘=1−𝑄Δ⊕𝐷₀⁻¹𝐷_𝑘+1−𝛾𝐹_𝑘𝑄₃×^∞₀𝑒^Δ𝑢⊗𝐼_𝑊^𝑢₀𝐼_𝑁+1⊗𝑃𝑖−𝑘,𝑦𝑒^{−𝑠(𝑢−𝑦)}(𝑢)+𝑑𝑦𝑑𝐵_𝑖𝑗=1𝝅_𝑗𝑄₁+𝛾𝑄₂^∞₀𝑒^Δ𝑢⊗𝐼_𝑊^𝑢₀𝐼_𝑁+1⊗𝑃𝑖−𝑗,𝑦𝑒^{−𝑠(𝑢−𝑦)}(𝑢)+𝑑𝑦𝑑𝐵1−𝛾^𝑖−𝑗_𝑘=0𝐹_𝑘𝑄₃^∞₀𝑒^Δ𝑢⊗𝑃𝑖−𝑘−𝑗,𝑦𝑒^{−𝑠(𝑢−𝑦)}𝑑𝑦𝑑𝐵(𝑢),𝐯(𝑖,2,𝑠)=𝜏⁻¹1−𝛾^𝑖−1_𝑗=0𝝅_𝑗^∞₀𝑑𝐹(𝑢)⊗𝐼_𝑊𝑄₃×^𝑢₀𝑒^{−𝑠(𝑢−𝑦)}𝐼_𝑁+1⊗𝑃𝑖−𝑗−1,𝑦𝑑𝑦𝐵^∗(𝑠)⊗𝐼_𝑊,𝑖>0. (4.2)(4.3)(4.4)(4.5)

Proof. Let 𝜅_𝑗^{(𝑖,𝑚,𝑥,𝑡)}(𝑟,𝜈;𝑟^′,𝜈^′) denote the conditional probability that, given time 0 is an instant of the service completion at Station 1 and the embedded Markov chain 𝜉_𝑛 is in the state (𝑗,𝑟,𝜈) at that time, the next service completion epoch at Station 1 occurs later than 𝑡, the discrete components of the process 𝜒_𝑡 take values (𝑖,𝑚,𝑟′,𝜈′) at time 𝑡 and the continuous-time component 𝑣_𝑡<𝑥.
Let us arrange the probabilities 𝜅_𝑗^{(𝑖,𝑚,𝑥,𝑡)}(𝑟,𝜈;𝑟^′,𝜈^′), for fixed values 𝑖,𝑗,𝑚, according to the lexicographic order of the states (𝑟,𝜈;𝑟′,𝜈′) and form the square matrices 𝐾_𝑗(𝑖,𝑚,𝑥,𝑡)=𝜅_𝑗^{(𝑖,𝑚,𝑥,𝑡)}𝑟,𝜈;𝑟^′,𝜈^′_{𝜈,𝜈′=0,𝑊;𝑟,𝑟′=0,𝑁}. (4.6)
Then, using the ergodic theorem for semi-regenerative processes, (see [17, Theorem  6.12]), the probability vectors 𝐕(𝑖,𝑚,𝑥) can be related to the stationary distribution 𝝅_𝑗,𝑗≥0, of the embedded Markov chain 𝜉_𝑛,𝑛≥1, by 𝐕(𝑖,𝑚,𝑥)=𝜏^∞−1_𝑗=0𝝅_𝑗^∞₀𝐾_𝑗(𝑖,𝑚,𝑥,𝑡)𝑑𝑡,𝑖≥0,𝑚=0,2. (4.7) The corresponding vector Laplace-Stieltjes transforms 𝐯(𝑖,𝑚,𝑠) are defined by 𝐯(𝑖,𝑚,𝑠)=𝜏^∞−1_𝑗=0𝝅_𝑗^∞₀𝐾^∗_𝑗(𝑖,𝑚,𝑠,𝑡)𝑑𝑡,𝑖≥0,𝑚=0,2, (4.8) where 𝐾^∗_𝑗(𝑖,𝑚,𝑠,𝑡)=∫^∞₀𝑒^−𝑠𝑥𝑑𝐾_𝑗(𝑖,𝑚,𝑥,𝑡).
From (4.8), formulas (4.2) follow immediately when we note that 𝐾^∗_𝑗(𝑖,𝑚,𝑠,𝑡)=0 for the range of arguments {𝑗≥0,𝑖=0,𝑚=1,2} and {𝑗≥0,𝑖>0,𝑚=0}.
Let 𝑚=1. The lengthy but straightforward expressions for the matrices 𝐾^∗_𝑗(𝑖,1,𝑠,𝑡), 𝑖>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 𝐯(𝑖,1,𝑠),𝑖>0. Similar calculations yield (4.3) and (4.5).

Further, we study the generalized service time distribution at Station 1.

Let 𝐵(𝑥) be the matrix distribution function of generalized service time. More specifically, let 𝐵(𝑥)=(𝐵(𝑥)_𝑟,𝑟′)_{𝑟,𝑟′=0,𝑁}, where 𝐵(𝑥)_𝑟,𝑟′=𝑃{𝑡_𝑛+1−𝑡_𝑛<𝑥,𝑟_𝑛+1=𝑟′∣𝑟_𝑛=𝑟,𝑖_𝑛≠0}. Denote ℬ(𝑠)=∫^∞₀𝑒^−𝑠𝑡𝑑𝐵(𝑡).

Lemma 4.2. The matrix Laplace-Stieltjes transform of the generalized service time distribution at Station 1 is calculated as ℬ(𝑠)=𝑄₁+𝛾𝑄₂+1−𝛾𝐹^∗(𝑠)𝑄₃𝐵^∗(𝑠), (4.9) where 𝐹^∗(𝑠)=∫^∞₀𝑒^−𝑠𝑡𝑑𝐹(𝑡).

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 (𝑄₁+𝛾𝑄₂)𝐵^∗(𝑠). 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−𝛾)∫^∞₀𝑒^−𝑠𝑡𝑑𝐹(𝑡)𝑄₃𝐵^∗(𝑠). The stated result (4.9) follows immediately.

Now we are ready to derive the equation for the vector Laplace-Stieltjes transform 𝐯₁(𝑠) of the distribution of the virtual sojourn time at Station 1. Let 𝑣₁(𝑟,𝜈,𝑥) be the probability that, at an arbitrary epoch, the BMAP is in state 𝜈, the virtual sojourn time at Station 1 is less than 𝑥, and the number of busy servers at Station 2 just before the end of the virtual sojourn time is 𝑟. Then 𝐯₁(𝑠) is defined as a vector of Laplace-Stieltjes transforms 𝑣₁(𝑟,𝜈,𝑠)=∫^∞₀𝑒^−𝑠𝑥𝑑𝑣₁(𝑟,𝜈,𝑥) written in lexicographic order.

Theorem 4.3. The vector Laplace-Stieltjes transform 𝐯₁(𝑠) satisfies the equation 𝐯₁(𝑠)𝐴(𝑠)=𝝅₀Φ(𝑠), (4.10) where 𝐴(𝑠)=𝑠𝐼+^∞_𝑟=0ℬ^𝑟(𝑠)⊗𝐷_𝑟,Φ(𝑠)=𝜏⁻¹𝑄Δ⊕𝐷₀⁻¹Δ⊗𝐼_𝑊−𝑠𝐼𝐵^∗(𝑠)⊗𝐼_𝑊. (4.11)

Proof. As mentioned above, the virtual sojourn time at Station 1 consists of the residual time from an arbitrary time 𝑡 to the next service completion epoch, the generalized service times of customers that await for a service at time 𝑡, 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 transform 𝐯₁(𝑠) as follows:𝐯₁(𝑠)=𝐯(0,0,𝑠)+^∞_𝑖=1𝐯(𝑖,1,𝑠)ℬ^𝑖(𝑠)⊗𝐼_𝑊+^∞_𝑖=1𝐯(𝑖,2,𝑠)ℬ^𝑖−1(𝑠)⊗𝐼_𝑊. (4.12)
Further, we multiply (4.12) by the matrix 𝑠𝐼+∑^∞_𝑟=0ℬ^𝑟(𝑠)⊗𝐷_𝑟, and, after some laborious calculations, get𝐯₁(𝑠)𝑠𝐼+^∞_𝑟=0ℬ^𝑟(𝑠)⊗𝐷_𝑟=𝜏⁻¹𝝅₀^∞_𝑖=0𝐶_𝑖+^𝑖+1_𝑗=1𝝅_𝑗𝑌_{𝑖−𝑗+1}ℬ^𝑖+1(𝑠)⊗𝐼_𝑊+𝝅₀𝑄−Δ⊕𝐷₀⁻¹𝑠𝐼+𝐷₀+𝐼𝐵^∗(𝑠)⊗𝐼_𝑊−^∞_𝑗=0𝝅_𝑗ℬ^𝑗+1(𝑠)⊗𝐼_𝑊. (4.13) Multiplying the balance equations for stationary probability vectors 𝝅_𝑖 of the form 𝝅_𝑖=𝝅₀𝐶_𝑖+^𝑖+1_𝑙=1𝝅_𝑙𝑌_{𝑖−𝑙+1},𝑖≥0, (4.14) by ℬ^𝑖+1(𝑠)⊗𝐼_𝑊 and summing over 𝑖 we obtain ^∞_𝑖=0𝝅_𝑖ℬ^𝑖+1(𝑠)⊗𝐼_𝑊=𝝅₀^∞_𝑖=0𝐶_𝑖ℬ^𝑖+1(𝑠)⊗𝐼_𝑊+^∞_𝑖=0^𝑖+1_𝑗=1𝝅_𝑗𝑌_{𝑖−𝑗+1}ℬ^𝑖+1(𝑠)⊗𝐼_𝑊. (4.15) Using (4.15) to simplify equation (4.13), we obtain (4.10).

4.1.2. The Sojourn Time at Station 2

Let 𝐯₂(𝑠) be the column vector of the Laplace-Stieltjes transforms of the conditional sojourn time distributions at Station 2. The 𝑟th 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 𝑟 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 𝐯₂(𝑠)=𝑄₄𝐹^∗(𝑠)+𝐼𝐼̂+𝛾𝑄₂+1−𝛾𝐹^∗(𝑠)diag𝑓_𝑟,0(𝑠),𝑟=𝑁,𝑁−1,…,0𝑄₃𝐞. (4.16)

Proof. The sojourn time of a customer who requires 𝑚 servers at Station 2 consists of:(i)the service time of the customer when at least 𝑚 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 type 𝑚 is performed by 𝑚 servers independently of each other and finishes when all 𝑚 servers complete the service. The distribution of this service time is defined by the Laplace-Stieltjes transform 𝑓_𝑚,0(𝑠),𝑚=1,𝑁.
Taking into account this fact together with (i)–(iii) and using the matrix notation, we obtain expression (4.16) for the vector 𝐯₂(𝑠). 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 𝑣(𝑠)=𝐯₁(𝑠)𝐼_𝑁+1⊗𝐞_𝑊𝐯₂(𝑠), (4.17) where the vectors 𝐯₁(𝑠) and 𝐯₂(𝑠), 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

Let 𝑣₁^(𝑎)(𝑠) and 𝑣^(𝑎)(𝑠) 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: 𝑣₁^(𝑎)(𝑠)=𝜆⁻¹𝐯₁(𝑠)^∞_𝑘=0ℬ^𝑘(𝑠)(ℬ(𝑠)−𝐼)⁻¹⊗𝐷_𝑘𝐞. (4.18)

Proof. The actual sojourn time at Station 1 of an arbitrary-tagged customer, who arrived in a group of size 𝑘 and placed at the 𝑗th 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 𝑗−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 the 𝑘-size group that contains the tagged customer is evidently given by the vector 𝐯₁(𝑠)(𝐼_𝑁+1⊗𝑘𝐷_𝑘𝐞/𝜆).
Assuming that an arbitrary customer arriving in a group of size 𝑘 is placed on the 𝑗th position with probability 1/𝑘 and using the law of total probability, we immediately obtain the following expression: 𝑣₁^(𝑎)(𝑠)=^∞_𝑘=1𝐯₁(𝑠)𝐼_𝑁+1⊗𝑘𝐷_𝑘𝐞𝜆^𝑘_𝑗=11𝑘ℬ^𝑗−1(𝑠)𝐞. (4.19) 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: 𝑣^(𝑎)(𝑠)=𝜆⁻¹𝐯₁(𝑠)^∞_𝑘=0ℬ^𝑘(𝑠)(ℬ(𝑠)−𝐼)⁻¹⊗𝐷_𝑘𝐼_𝑁+1⊗𝐞_𝑊𝐯₂(𝑠). (4.20)

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 𝑠=0. This requires the calculation of the derivatives of 𝐯₁(𝑠) at the point 𝑠=0. However, the matrix 𝐴(𝑠) in (4.10) is singular at the point 𝑠=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 notation 𝐯₁^(𝑚)(𝑠) for the 𝑚th derivative of the vector 𝐯₁(𝑠),𝑚≥1,and set 𝐯₁⁽⁰⁾(𝑠)=𝐯₁(𝑠). Similar notations will be used for other functions of 𝑠.

Theorem 4.8. Let ∫^∞₀𝑡^𝑚𝑑𝐵(𝑡)<∞,𝑚=1,𝑀+1, where 𝑀 is an arbitrary positive integer. Then the vectors 𝐯₁^(𝑚)(0),𝑚=1,𝑀, are computed recursively by 𝐯₁^(𝑚)(0)=𝝅₀Φ^(𝑚)(0)−^𝑚−1_𝑙=0𝑚𝑙𝐯₁^(𝑙)(0)𝐴^{(𝑚−𝑙)}(0)𝐼̃+1𝝅𝑚+1⁰Φ^(𝑚+1)(0)−^𝑚−1_𝑙=0𝑙𝐯𝑚+1_1(𝑙)(0)𝐴^{(𝑚+1−𝑙)}(0)𝐞𝐞̂𝐴⁻¹, (4.21) with the initial condition 𝐯₁⁽⁰⁾(0)=𝐯₁(0)=𝜏⁻¹𝝅₀𝑄Δ⊕𝐷₀⁻¹Δ𝐵^∗(0)⊗𝐼_𝑊𝐼̃+𝑝_idle𝐞̂𝐴⁻¹, (4.22) where 𝐴=𝐴(0)𝐼̃+𝐴′(0)𝐞𝐞̂.

Proof of the theorem is presented in Appendix B. In what follows we assume that ∫^∞₀𝑡^𝑘𝑑𝐵(𝑡)<∞,𝑘=1,2.

Corollary 4.9. The mean virtual sojourn time at Station 1 is given by 𝑣₁=𝜏⁻¹𝝅₀𝑄Δ⊕𝐷₀⁻¹𝐵^∗(0)−Δ𝐵^∗′(0)⊗𝐼_𝑊+𝐯₁(0)𝐴^′(0𝐼̃+𝑝_idle)𝑏₁+12𝐯₁(0)𝐴^′′(0)𝐞𝐞̂𝐴⁻¹𝐞, (4.23) where the vector 𝐯₁(0) is defined by formula (4.22).

Proof follows from the relation 𝑣₁=−𝐯^′₁(0)𝐞 and formula (4.21).

Theorem 4.10. The mean virtual sojourn time in the system is given by 𝑣=𝑣₁+𝐯₁(0)𝐼_𝑁+1⊗𝐞_𝑊𝐯₂, (4.24) where 𝐯₁(0) and 𝑣₁ are given in (4.22) and (4.23), respectively, and 𝐯₂ is a vector of conditional means of the sojourn time at Station 2, 𝐯₂=−𝑄₄𝐹⁽¹⁾𝐞̂^𝑇+1−𝛾𝐹⁽¹⁾−𝐸diag^𝑟_𝑙=1𝑙𝜇⁻¹,𝑟=𝑁,𝑁−1,…,0𝑄₃𝐞, (4.25) where 𝐹⁽¹⁾ is given by formula (C.3) below.

Proof. To calculate the value 𝑣 we differentiate the expression in (4.17). Setting 𝑠=0 and replacing the sign, we have that 𝑣=−𝐯^′₁(0)𝐼_𝑁+1⊗𝐞_𝑊𝐯₂(0)−𝐯₁(0)𝐼_𝑁+1⊗𝐞_𝑊𝐯^′₂(0). (4.26) Putting 𝑠=0 in (4.16) we get 𝐯₂(0)=[𝑄₄(𝐸+𝐼)𝐼̂+𝛾𝑄₂+(1−𝛾)𝐸𝑄₃]𝐞=𝐞. This implies that the first term in the right-hand side of (4.26) is equal to −𝐯^′₁(0)𝐞=𝑣₁. Using the relation 𝐯₂=−𝐯^′₂(0) and differentiating (4.16) at the point 𝑠=0, we readily verify that 𝐯₂ has form (4.25). This completes the proof.

Theorem 4.11. The mean actual sojourn time at Station 1 is given by 𝑣₁^(𝑎)=−𝜆⁻¹𝐯^′₁(0)𝐞⊗^∞_𝑘=1𝑘𝐷_𝑘𝐞+𝐯₁(0)^∞_𝑘=1𝐼_𝑁+1⊗𝐷_𝑘𝐞^𝑘−1_𝑛=1^𝑛−1_𝑙=0ℬ^𝑙(0)ℬ^′(0)𝐞. (4.27)

Proof of the theorem follows from the relation 𝑣₁^(𝑎)=−𝑑𝑣₁^(𝑎)(𝑠)/𝑑𝑠|_𝑠=0 and formula (4.19).

Corollary 4.12. The mean actual sojourn time in the system is given by 𝑣_𝑎=𝑣₁^(𝑎)+𝜆⁻¹𝐯₁(0)^∞_𝑘=1𝐼_𝑁+1⊗𝐷_𝑘𝐞^𝑘−1_𝑛=0ℬ^𝑛(0)𝐯₂. (4.28)

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 time 𝑏₁=0.1, but different values for the coefficient of variation, 𝑐_var. The first process is coded as 𝐷 and corresponds to the deterministic service time distribution. The second process is coded as 𝑀 and corresponds to the exponential service time. The third and the fourth service processes are coded as 𝐻𝑀₂⁽¹⁾, 𝐻𝑀₂⁽²⁾ 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 𝐻𝑀₂⁽¹⁾ and (0.98,0.02) and the intensities 10000,0.2 in the case of 𝐻𝑀₂⁽²⁾. The coefficients of variation of processes 𝐷, 𝑀, 𝐻𝑀₂⁽¹⁾, and 𝐻𝑀₂⁽²⁾ are, respectively, equal to 0,1,5,9.95.

The input process is defined by the matrices𝐷₀=.−15.73270.60620.59240.5178−2.28970.46790.59710.5653−1.9597,𝐷=14.15020.30210.08180.10711.0320.16460.08580.19790.5136 (5.1) The matrices 𝐷_𝑘,𝑘=1,5, are calculated as follows: 𝐷_𝑘=𝐷ℎ^𝑘−1(1−ℎ)/(1−ℎ⁵), where ℎ=0.8. Then we normalize the matrices 𝐷_𝑘,𝑘=0,5, so as to get the arrival rate 𝜆=1. This BMAP has coefficient of correlation 𝑐_cor=0.2. The other parameters of the system are as follows: 𝑁=5,𝜇=0.8,𝛾=0.5,𝑞₀=0.1,𝑞₁=𝑞₂=0.3,𝑞₃=𝑞₄=𝑞₅=0.1.

We vary the mean service time 𝑏₁ 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 𝑏₁ and 𝑐_var. The values of 𝜌 that exceed 1 are printed in bold face. The tandem queueing system is not stable for these values.

Figures 1 and 2 show the dependence of the main performance measures of the system on the value of the mean service time 𝑏₁ 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.

(a)

(b)

Figure 1 

The mean virtual and actual sojourn time as functions of the mean service time for different service time distributions.

(a)

(b)

Figure 2 

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.

Experiment 2. Here we solve numerically the following optimization problem. Find an optimal choice for the number 𝑁 of servers at Station 2 that will minimize the expected total cost per unit of time: 𝐽=𝐽(𝑁)=𝑎𝑁+𝑐₁𝜆𝑃_loss+𝑐₂𝑣^(𝑎), (5.2) where 𝑎 is the cost of utilization per unit time of a server at Station 2 (maintenance cost), 𝑐₁ is the cost of a customer leaving the system after a service at Station 1 due to lack of required servers, and 𝑐₂ is the cost per unit of time of holding (the sojourn time) an arbitrary customer in the system (holding cost).

Using the MAP characterized by the matrices𝐷₀=−6.745385.45412×10⁻⁶5.45412×10⁻⁶,−0.219455,𝐷=6.7006850.0446950.1224270.097023 (5.3) we construct the BMAP with the matrices 𝐷_𝑘,𝑘=0,5, similar to Experiment 1 and normalize so as to have 𝜆=3. The BMAP has the coefficient of correlation 𝑐_cor=0.2 and the coefficient of variation 𝑐_var=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 𝐪 are: 𝑞₀=0.1,𝑞₁=0.9,𝑞_𝑚=0,𝑚=2,𝑁. The various costs are taken as follows: 𝑎=5,𝑐₁=50,𝑐₂=3.

The objective function, 𝐽, as a function of the number 𝑁 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.

Figure 3 

The objective function as a function of the number of servers in Station 2 for different service rates.

The optimal values 𝐽^∗ 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 𝑁^∗ increases from 3 to 6. The relative gain of the optimal configuration in comparison to a system with an arbitrary number 𝑁 of the servers at Station 2 is defined as 𝑅_rel(𝑁)=((𝐽(𝑁)−𝐽^∗)/𝐽^∗)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 𝐽^∗ is 20.6657 and the optimal number of servers 𝑁^∗=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 𝑁^∗=3 instead of 2 servers and maximal relative gain is more than 161% if we use 𝑁^∗=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 correlation 𝑐_cor=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 𝐷₀ and 𝐷₁=𝐷, from which the matrices 𝐷_𝑘,𝑘=0,5, are defined in the same way as in Experiment 1.

The BMAP having the coefficient of correlation 𝑐_cor=0.1 is characterized by the matrices𝐷₀=.−13.33460.58860.61730.6927−2.44660.42290.68230.4144−1.6354,𝐷=11.54690.36310.21870.38420.86590.08090.28520.04250.2111 (5.4) The BMAP having the coefficient of correlation 𝑐_cor=0.3 is defined by the matrices𝐷₀=.−25.53980.39330.36120.1452−2.23220.20000.29600.3874−1.7526,𝐷=24.24210.46690.07630.03411.66680.18610.00900.25550.8047 (5.5) All the three defined BMAPs have the coefficient of variation 𝑐_var=2.

In addition, we consider the BMAP which is a group Poisson process. It has the same mean arrival rate as three other BMAPs, coefficient of correlation 𝑐_cor=0 and the coefficient of variation 𝑐_var=1.

The service time at Station 1 has the Erlangian distribution of order 3 with parameter 20. The mean service time 𝑏₁=3/20 and the squared coefficient of variation 𝑐²_var=1/3.

The number of servers at Station 2 𝑁=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 𝐪, which defines the type of the customer, are 𝑞₀=0.1,𝑞₁=𝑞₂=0.3, and 𝑞₃=𝑞₄=𝑞₅=0.1.

Figures 4 and 5 illustrate the dependence of the mean virtual sojourn time 𝑣, the mean actual sojourn time 𝑣_𝑎, the loss probability 𝑃_loss, and the mean number 𝑁_busy 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.

(a)

(b)

Figure 4 

The mean virtual and actual sojourn time as functions of the system load for the BMAPs with different correlation.

(a)

(b)

Figure 5 

The loss probability and the mean number of busy servers at Station 2 as functions of the system load for the BMAPs with different correlation.

Figures 4 and 5 confirm the fact that values of 𝑣, 𝑣_𝑎, 𝑃_loss, and 𝑁_busy 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, the BMAP/𝐺/1→•/𝑀/𝑁/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.

Appendices

A. Expressions for the Matrices 𝐾^∗_𝑗(𝑖,1,𝑠,𝑡),𝑗>0,𝑖>0

𝐾^∗₀(𝑖,1,𝑠,𝑡)=𝑄₁+𝛾𝑄₂^𝑡₀𝐼_𝑁+1⊗𝑒^{𝐷0𝑥𝑖}_𝑘=1𝐷_𝑘(×𝑑𝑥𝑃𝑖−𝑘,𝑡−𝑥)_∞0𝑒^−𝑠𝑢𝑒^{Δ(𝑡+𝑢)}⊗𝐼_𝑊()+×𝑑𝐵𝑡−𝑥+𝑢1−𝛾_𝑡0𝑑𝐹(𝑥)⊗^𝑖_𝑘=1𝑃(𝑘,𝑥)𝑃(𝑖−𝑘,𝑡−𝑥)𝑄₃×^∞₀𝑒^−𝑠𝑢𝑒^{Δ(𝑡−𝑥+𝑢)}⊗𝐼_𝑊𝑑𝐵(𝑡−𝑥+𝑢)+^𝑡₀^𝑦₀𝑑𝐹(𝑥)⊗𝑒^{𝐷0𝑦𝑖}_𝑘=1𝐷_𝑘𝑑𝑦𝑃𝑖−𝑘,𝑡−𝑦𝑄₃×^∞₀𝑒^−𝑠𝑢𝑒^{Δ(𝑡−𝑥+𝑢)}⊗𝐼_𝑊𝐾𝑑𝐵𝑡−𝑦+𝑢,𝑖>0,_∗𝑗(𝑖,1,𝑠,𝑡)=𝑄₁+𝛾𝑄₂𝐼_𝑁+1⊗𝑃𝑖−𝑗,𝑡^∞₀𝑒^−𝑠𝑢𝑒^{Δ(𝑡+𝑢)}⊗𝐼_𝑊()+𝑑𝐵𝑡+𝑢1−𝛾^𝑡₀𝑑𝐹(𝑥)⊗^𝑖−𝑗_𝑘=0𝑃(𝑘,𝑥)𝑃𝑖−𝑘−𝑗,𝑡−𝑥𝑄₃×^∞₀𝑒^−𝑠𝑢𝑒^{Δ(𝑡−𝑥+𝑢)}⊗𝐼_𝑊𝑑𝐵(𝑡−𝑥+𝑢). (A.1)

B. Proof of Theorem 4.8

Proof. Successively differentiating the expression in (4.10), we get𝐯₁^(𝑚)(0)𝐴(0)=𝝅₀Φ^(𝑚)(0)−^𝑚−1_𝑙=0𝑚𝑙𝐯₁^(𝑙)(0)𝐴^{(𝑚−𝑙)}(0),𝑚≥0. (B.1) It follows from (4.11) that 𝐴(0)=∑^∞_𝑟=0ℬ^𝑟(0)⊗𝐷_𝑟, where ℬ(0) is an irreducible stochastic matrix. This implies that 𝐴(0) is an irreducible infinitesimal generator, and hence 𝐴(0) is a singular matrix. Thus, it is not possible to develop a recursive scheme for computing the vectors 𝐯₁^(𝑚)(0),𝑚≥0, directly from (B.1). We will now modify the system (B.1) to get the system with a nonsingular matrix. To this end, we postmultiply the expression for 𝑚+1 in (B.1) on both sides with 𝐞. Taking into account that 𝐴(0)𝐞=𝟎^𝑇, we get𝐯₁^(𝑚)𝐴^′(0)1𝐞=𝝅𝑚+1⁰Φ^(𝑚+1)(0)−^𝑚−1_𝑙=0𝑙𝐯𝑚+1_1(𝑙)(0)𝐴^{(𝑚+1−𝑙)}(0)𝐞. (B.2) It can be shown that the right-hand side of (B.2) is not equal to zero. It is positive if 𝑚=2𝑘 and it is negative if 𝑚=2𝑘+1,𝑘≥0. Thus, replacing one of the equations in the system (B.1) (without loss of generality, we replace the first equation) with equation (B.2), we get the following (inhomogeneous) system of linear algebraic equations for the entries of the vector 𝐯₁^(𝑚)(0)∶𝐯₁^(𝑚)(0)𝐴=𝝅₀Φ^(𝑚)(0)−^𝑚−1_𝑙=0𝑚𝑙𝐯₁^(𝑙)(0)𝐴^{(𝑚−𝑙)}(0)𝐼̃+1𝝅𝑚+1⁰Φ^(𝑚+1)(0)−^𝑚−1_𝑙=0𝑙𝐯𝑚+1_1(𝑙)(0)𝐴^{(𝑚+1−𝑙)}(0)𝐞𝐞̂,𝑚≥0. (B.3) The above system has the unique solution if 𝐴 is non-singular. We prove this by showing that det𝐴≠0.

Let us calculate det𝐴 asdet𝐴=∇𝐴^′(0)𝐞, (B.4) where ∇ is a vector of algebraic cofactors of the first column of the matrix 𝐴(0). Since 𝐴(0) is irreducible, the vector ∇ is proportional to any solution of the system𝐱𝐴(0)=𝟎, (B.5) that is,∇=𝑐𝐱, (B.6) where the scalar 𝑐≠0.

Let the vector 𝝑 be the unique solution to the system𝝑ℬ(0)=𝜗,𝝑𝐞=1. (B.7) Then, by the direct substitution it can be verified that the vector 𝐱=𝜗⊗𝜽 is a solution to the system (B.5).

From (B.6), ∇=𝑐(𝝑⊗𝜽). Substituting the vector ∇ into (B.4), we obtaindet𝐴=𝑐(𝝑⊗𝜽)𝐴^′(0)𝐞=𝑐(𝝑⊗𝜽)𝑠𝐼+^∞_𝑟=0ℬ^𝑟(𝑠)⊗𝐷_𝑟^′_𝑠=0𝐞=𝑐+𝑐^∞_𝑟=1𝑟𝝑ℬ^′(0)𝐞⊗𝜽𝐷_𝑟𝐞=𝑐1+𝜆𝝑ℬ^′(0)𝐞. (B.8) In further evaluation of det𝐴, we use the ergodicity condition given in Theorem 3.3.

Setting 𝑠=0 in (4.9) and noting that 𝑄₁+𝛾𝑄₂+(1−𝛾)𝐹^∗(0)𝑄₃=𝑄, we see that ℬ(0)=𝑄𝐵^∗(0) and that the vector 𝜗 defined by (B.7) is the unique solution of system (3.8). It can be easily verified that𝝑ℬ^′(0)𝐞=−𝑏₁+1−𝛾𝜗^∞₀𝑡𝑑𝐹(𝑡)𝑄₃𝐞. (B.9) Multiplying (B.9) by 𝜆 and comparing the obtained equation with the equality in (3.12) we see that𝜆𝝑ℬ^′(0)𝐞=−𝜌. (B.10) It follows from (B.8) and (B.10) that det𝐴=𝑐(1−𝜌). Since the stability condition implies that 𝜌<1 and since 𝑐≠0, we have det𝐴≠0. This competes the proof of the theorem.

C. Formulas for Calculation of the Mean Virtual Sojourn Time

The expressions (4.21)-(4.22) for the mean virtual sojourn time involve the matrices 𝐴,𝐴′(0) and the vector 𝐴^′′(0)𝐞 for which we derive explicit expressions below.𝐴=^∞_𝑟=0ℬ^𝑟(0)⊗𝐷_𝑟𝐼̃+𝐼+^∞_𝑟=1^𝑟−1_𝑛=0ℬ^𝑛(0)ℬ^′(0)⊗𝐷_𝑟𝐞̂𝐴𝐞,_′(0)=𝐼+^∞_𝑟=1^𝑟−1_𝑛=0ℬ^𝑛(0)ℬ^′(0)ℬ^{𝑟−𝑛−1}(0)⊗𝐷_𝑟,𝐴^′′(0)𝐞=^∞_𝑟=1^𝑟−1_𝑛=12^𝑛−1_𝑙=0ℬ^𝑙(0)ℬ^′(0)ℬ^{𝑛−𝑙−1}(0)ℬ^′(0)+ℬ^𝑛(0)ℬ^′′(0)⊗𝐷_𝑟𝐞, (C.1) where ℬ(0)=𝑄𝐵^∗(0),ℬ^′(0)=𝑄𝐵^∗(1)+1−𝛾𝐹⁽¹⁾𝑄₃𝐵^∗(0),ℬ^′′(0)=𝑄𝐵^∗(2)−21−𝛾𝐹⁽¹⁾𝑄₃𝐵^∗(1)+1−𝛾𝐹⁽²⁾𝑄₃𝐵^∗(0),𝐹^(𝑚)=(−1)^𝑚∞₀𝑡^𝑚𝑑𝐹(𝑡)=𝑓^(𝑚)_{𝑟,𝑟′𝑟,𝑟′=0,𝑁},𝑚=1,2,𝑓⁽¹⁾_𝑟,𝑟′=−0,𝑟≤𝑟′,_{𝑟𝑙=𝑟′+1}𝑙𝜇⁻¹,𝑟>𝑟′,𝑓⁽²⁾_𝑟,𝑟′=20,𝑟≤𝑟′,_{𝑟𝑙=𝑟′+1}(−1)^{𝑙−𝑟′+1}𝑟−𝑟′𝑙−𝑟′𝑙𝜇⁻²𝐵,𝑟>𝑟′,_∗(𝑚)=(−1)^𝑚∞₀𝑡^𝑚𝑒^Δ𝑡𝑑𝐵(𝑡)=𝛽^∗(𝑚)_{𝑟,𝑟′𝑟,𝑟′=0,𝑁},𝑚=𝛽0,2,_{∗(𝑚)𝑟,𝑟′}=𝑟0,𝑟<𝑟′,𝑟′^{𝑟−𝑟′}_𝑖=0(−1)^𝑖𝑟−𝑟^′𝑖𝛽^(𝑚)𝜇𝑟^′+𝑖,𝑟≥𝑟′. (C.2)(C.3)(C.4)(C.5)

Remark C.1. Formulas for 𝐴,𝐴^′(0),𝐴^′′(0) contain infinite sums. However, the calculations of these should not create any difficulty as for overwhelming majority of interesting and useful queueing models; the parameter matrices, 𝐷_𝑘, of the BMAP process are equal to zero for 𝑘 greater than some threshold, say, 𝐾. Thus, all sums become finite. In alternative case, some analytical formula for computing the infinite sequence of matrices 𝐷_𝑘,𝑘≥1, should be given. If this sequence is generated such as 𝐷_𝑘=𝐷(1−𝜎)𝜎^𝑘−1,𝑘≥1, where 𝐷 and 0<𝜎<1 are given parameters, the infinite sums can be computed explicitly.

Acknowledgments

The authors are thankful to Professor S. R. Chakravarthy (Kettering University, USA) for valuable remarks and suggestions and kind help in polishing the text of this paper. This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (Grant no. 2010-0003269).