ALGORITHMIC ANALYSIS OF THE MAXIMUM LEVEL LENGTH IN GENERAL-BLOCK TWO-DIMENSIONAL MARKOV PROCESSES

Two-dimensional continuous-time Markov chains (CTMCs) are useful tools for studying stochastic models such as queueing, inventory, and production systems. Of particular interest in this paper is the distribution of the maximal level visited in a busy period because this descriptor provides an excellent measure of the system congestion. We present an algorithmic analysis for the computation of its distribution which is valid for Markov chains with general-block structure. For a multiserver batch arrival queue with retrials and negative arrivals, we exploit the underlying internal block structure and present numerical examples that reveal some interesting facts of the system.


Introduction
It is well-known that two-dimensional Markov processes often provide an adequate initial starting point for describing the behavior of many stochastic models.The detailed analysis of particular structures has been the subject matter of many books and papers [8,12,13].Specifically, these address the importance of Markov chains whose transition probability matrices (discrete case) or infinitesimal generators (continuous case) are of quasi-birth-and-death (QBD), M/G/1, or GI/M/1 type.For such structures, the matrixanalytic methods provide implementable and numerically stable solutions for the distribution of the system state and other descriptors of the performance quality.
In this paper, we develop an algorithmic analysis of the maximum level visited by a two-dimensional CTMC during a busy period.Unlike the study of the stationary system state, the analysis of extreme values can be done for general-block structures.So, specific well-suited block forms are not assumed here.The analysis of extreme values often provides a performance measure of interest connected with the system congestion.In a queueing model, for example, moderate values of the queue indicate that the system is evolving smoothly, whereas large values indicate the need for one or more actions such as increasing the service capacity and rescheduling of work, on the part of the system manager.Extreme values of Markov processes can be investigated using different approaches.For example, in [11] the objective is to determine the distribution of the maximum queue length during a busy period, and in [14] an asymptotic approach is proposed.In this paper, we will adopt the former approach.
In addition to the assumption of dealing with a well-posed matrix structure, in the framework of the two-dimensional Markov chains, some typical requirements are the stationary regime and the system homogeneity.In the absence of these assumptions, the analysis becomes extremely intricate or, indeed, impossible.In the sequel, we will show how the maximum level visited in a busy period overcomes these restrictions and can be efficiently computed in a class of general-block CTMC with level-dependent structure operating in both stationary and nonstationary regimes.
On the other hand, in recent years, there has been an increasing interest in the investigation of retrial phenomenon in computer and telecommunication systems.The existing literature includes applications to collision-avoidance star protocols for local area networks (CASLANs) [7], circuit-switched networks with hybrid fiber-coax (HFC) architecture [6], fast reservation protocols for asynchronous transfer mode (ATM) networks [15], and cellular mobile telephone networks (CMTNs) [10,16].Thus, the interest of retrial queues in engineering is a motivating factor of this paper that focuses on applications and algorithmic implementation of such queueing models.
As related works, we mention [4] where a study for the maximum level visited is initiated for quasi-birth-and-death processes with applications to call centers.Extensions to other models involving Markovian arrival processes (MAPs) with multiple exponential servers [3] and Markovian arrivals with single-server and phase-type services [2] have recently been studied.
The paper is organized as follows.In Section 2, we describe the class of two-dimensional Markov chains under study.Then, in Section 3, we show how the computation of the maximum level visited by the process reduces to finding some absorption probabilities in an auxiliary CTMC with two absorbing states.To this end, we propose a "forwardelimination, backward-substitution" algorithm.In Section 4, we consider a multiserver retrial queue with batch arrivals which operates under the presence of a second flow of negative batch arrivals.We exploit the special nature of the involved blocks to refine the algorithm for the computation of the maximum number of customers in the retrial queue.Some numerical work for this queueing model is given in Section 5.More specifically, we illustrate the effect of the system parameters on the percentiles and modes of the maximum level visited.

The general-block two-dimensional Markov chain
We consider a CTMC {(N (t),C(t)); t ≥ 0} with state space given by S = ∪{(i, j) : i ≥ 0, 0 ≤ j ≤ L i }.In the state (i, j), the first coordinate i is called the level of the state.For use in the sequel, we will denote by 0 * the state (0,0).The other vector states are defined as follows: (2.1) J. R. Artalejo and S. R. Chakravarthy 3 We assume that the number of states in each level, L i , is finite.The infinitesimal generator of the CTMC, Q = [q (i, j)(m,n) ] (for later use, we define q (i, j) = −q (i, j)(i, j) ), is of the form We further assume that our CTMC is irreducible and regular.

Analysis of the maximal level visited in a busy period
We now define the busy period for the class of general-block two-dimensional Markov chains described in Section 2. In this context, the busy period is defined as the duration that starts when the process leaves state 0 * and ends at the first epoch thereafter when the process visits state 0 * again.
In queueing theory, the first transition after leaving state 0 * usually corresponds to an arrival of a customer.However, a concrete definition of a busy period depends on the specific stochastic model under study.For example, if customers arrive according to an MAP [5,9], then the first transition does not necessarily represent an arrival.In fact, stochastic systems modelled with MAP input often need a higher-dimensional representation.In these cases, the definition of busy period must be modified and the analysis in the sequel remains valid with appropriate modifications (see [2,3]).
Define N max as the maximum level visited by the Markov chain during a busy period.We first observe that the mass function of N max can be expressed as where S 0 = {(m, n) : q 0 * , (m,n) > 0}, that is, S 0 represents the set of states that the process can visit when a busy period starts.
For initial states (m,n) such that m > k, we trivially have P (m,n) (N max = k) = 0. We also note that P (m,n) (N max < k) corresponds to the probability of the event that starting from the state (m,n) the Markov chain {(N (t),C(t)); t ≥ 0} will hit state 0 * before hitting level k, for k ≥ 1.Thus, the computation of P(N max < k), for k ≥ 1, reduces to finding the probability of absorption to a particular state in a finite-state CTMC with two absorbing states, say, 0 * and k * , and whose infinitesimal generator is given below: where the coefficient matrices appearing in (3.2) are as seen in (2.2), and Let T(k) denote the part of the generator Q(k) that corresponds to the levels 0,1,..., k − 1.Since the original CTMC is irreducible, this set of states is transient and consequently T(k) is invertible.Let T 0 (k) denote a column vector obtained by stacking the column vectors A 0,0 * ,A 1,0 * ,..., and A k−1,0 * , and let D 0 (k) = (A 0 * ,0 ,A 0 * ,1 ,...,A 0 * ,k−1 ).Then, it is easy to verify that Of course, once the cumulative mass function P(N max < k) is computed, we automatically determine the mass function P(N max = k).We also note that when L 0 = 0, the blocks corresponding to level 0 in (3.2) must be deleted.
It should be pointed out that the analysis of a CTMC with a finite number K of levels needs a minor modification.Following the above arguments, we may determine Here and in the sequel the notation " " stands for the transpose of a matrix.The system T(k)x (k) = −T 0 (k) can be solved using any all-purpose algorithm.For the sake of completeness, we next formulate a "forward-elimination, backward-substitution" algorithm (Algorithm 3.1) .
It should be noted that matrices Δ j j , 0 ≤ j ≤ k − 1, are invertible.The proof relies on an inductive argument that uses the existence of the inverse of T(k), for any k ≥ 1.
In the case of particular block structures (i.e., QBD-, M/G/1-, or GI/M/1-type generators), the "forward-elimination, backward-substitution" algorithm can be appropriately reformulated.Some simplifications are also obtained by exploiting specific sparse structures of the blocks.Examples can be seen in [2,3].

Application to the M [X] /M/c retrial queue with negative arrivals
We consider a multiserver queueing system where the arrivals occur in batches according to a Poisson process with rate λ.The batch sizes are i.i.d.random variables with probability mass function given by {a k : Customers are served by c homogeneous exponential servers each serving at a rate μ.Arriving bath of customers finding any free server will get into service immediately and any customer not getting into service will enter into an orbit of infinite size.Customers in orbit will compete for service by searching for a free server.The search (retrial) times are assumed to be exponentially distributed with parameter θ.In addition, we also assume a second flow of negative arrivals occurring in batches according to a Poisson process with parameter δ.The effect of these arrivals is seen only when all c servers are busy with at least one customer waiting in the orbit.The batch sizes of negative arrivals are i.i.d.
Step 1 (forward elimination): Step 2 (backward substitution): random variables with probability mass function given by {d k : When a negative arrival occurs with a batch of size k at a time when the orbit size is i with all c servers busy, then r = min{k,i} customers are removed from the orbit.These r customers are considered lost.Putting a k = 1, the model reduces to the retrial queue with negative arrivals investigated in [1].
The infinitesimal generator of the M [X] /M/c retrial queue with negative arrivals has a general-block structure with the following blocks: (e j denotes the column vector of appropriate dimension with 1 in the jth position and 0 elsewhere.)
A i,0 , i ≥ 2, is a rectangular matrix of dimension (c + 1) × c whose elements are all zero except the (c + 1,c)th entry which is given by δ * i , is a square matrix of dimension c + 1 whose elements are all zero except the (c + 1,c + 1)th entry which is given by δ i− j .
A i, j , i ≥ 1, j ≥ i + 1, is a square matrix of dimension c + 1 whose elements are all zero except the last column which has entries given by {λ c+ j−i ,λ c+ j−i−1 ,...,λ j−i }.
A 1,0 is a rectangular matrix of dimension (c + 1) × c, and A i,i−1 , i ≥ 2, is a square matrix of dimension c + 1 which are given by Although the computation of the maximum number of customers in orbit during a busy period can be done for the case of nonstationary regime, this situation is quite unusual because the number of customers in the system grows unbounded.Hence, we will concentrate mainly on the stationary case.The CTMC {(N (t),C(t)); t ≥ 0} modelling the M [X] /M/c retrial queue with negative arrivals is positive recurrent if and only if λa < cμ + δd, ( where a and d denote, respectively, the expected batch size of positive and negative arrivals. J. R. Artalejo and S. R. Chakravarthy 7 To prove the necessity of (4.5), we partition the state space into diagonal subsets D(k)= {(i, j) : i + j ≤ k}.By equating the flow rate in and out of the subset D(k) and adding over k, we get where P i j denotes the stationary probability that the system is in state (i, j).
The first term on the right-hand side of (4.6) is bounded by cμ.On the other hand, by interchanging the order of summation in the second term of the right-hand side and taking into account that ∞ i=k P ic < 1, we easily conclude that λa < cμ + δd.The sufficiency can be proved by appealing to the classic Foster criterion for the drifts of the embedded Markov chain {Z n ; n ≥ 0} at the transition epochs.We may consider the test function f (i, j) = i + h j.In the case when i ≥ 0 and 0 ≤ j ≤ c − 1, the mean drifts are given by and, for i ≥ 1 and j = c, Let h be any number in the interval (max{(λa − δd)/cμ,0},1).Then, it is clear that we may find an appropriate ε > 0 such that γ(i, j) < −ε (except perhaps for a finite number of states).This proves that condition (4.5) is also sufficient for the positive recurrence.
In what follows, we will describe an efficient algorithm for the computation of the probability P(N max < k) that is of interest in this paper.Before proceeding with the description of the algorithm, we register a number of observations to ease the understanding of the equations proposed.
(i) To minimize the storage requirements of the number of auxiliary matrices needed, the equations are modified accordingly.For example, where we need matrices of the form (−A ii ) −1 , 1 ≤ i ≤ k − 1, we simply store (−A k−1,k−1 ) −1 and modify the corresponding equations appropriately.
(ii) The equations are displayed in such a way that they are ready for numerical implementation.At this point, we can use a variety of numerical methods such as Gaussian elimination, LU factorization and aggregate/disaggregate method.
Observing that formula (3.4) where can be written in terms of matrices of dimension at most equal to c + 1 that are well suited for numerical implementation.We need to present separately the cases c = 1 and c ≥ 2, and the cases k = 1 and k ≥ 2.
Here we compute the probability explicitly as In this case, the equations are given by and, for 1 Here the probability is calculated explicitly as where a = (a 1 ,...,a c ).
J. R. Artalejo and S. R. Chakravarthy 9 Case 4. c ≥ 1 and k ≥ 2. In this case, the equations are given by where (4.18)

Numerical examples
In order to evaluate the performance of the M [X] /M/c retrial queue with negative arrivals, some computational experiments are performed.To this end, we will consider the following three distributions for the batch size.(a) Geometric (GEO): (b) Deterministic (DET): (c) Conditional Poisson (POI): In the following α 1 , D 1 , and ν 1 (α 2 , D 2 , and ν 2 ) will indicate parameters associated to the batch of positive (negative) arrivals.The first set of numerical results (see Examples 5.1-5.3)shows the influence of the system parameters on the distribution of the maximum orbit length for a model with geometric batch sizes of positive and negative arrivals.More concretely, we will focus on two measures: the 99th percentile and the mode(s) of N max .In Example 5.4, we analyze the effect of the batch size distributions.Example 5.5 shows how the probability mass function of N max behaves as a function of the traffic intensity ρ = λa(cμ + δg) −1 .Finally, in Example 5.6, we consider the M [X] /M/c retrial queue (case δ = 0) and M/M/c retrial queue with negative arrivals (case a 1 = 1) which are examples, respectively, of structures of the M/G/1 and GI/M/1 types.
Example 5.1.The main purpose of this example is to see the influence of μ and θ on the two measures: the 99th percentile and the mode(s) of the random variable N max .In this example, we fix c = 5, λ = 1.5, α 1 = 0.5, δ = 1.75, α 2 = 0.3 and calculate the probability mass function of N max .μ is varied in such a way that ρ will take the values 0.2, 0.4, 0.6, and 0.8.In Table 5.1 we display the percentile and the mode(s) as ordered set of values.For example, the entry (10,0) will mean that the 99th percentile occurs at k = 10 and there is only mode that occurs at 0. The entry (45,0,28) indicates that the 99th percentile occurs at k = 45 and the two modes are at k = 0 and k = 28.
The following observations can be inferred from Table 5.1.
(i) For all values of μ, the 99th percentile appears to decrease as θ increases and seems to approach a constant value.This is to be expected since an increase in θ will result in more success for orbit customers to occupy a free server.
(ii) For all values of θ, the 99th percentile appears to increase as μ decreases.Again, this is as to be expected since a decrease in μ will result in seeing more customers in the orbit.
(iii) In all the cases, there appears to be a mode at 0 and depending on the values of μ and θ there seems to be multiple modes.In the case of multiple modes, it is interesting to see that the second mode appears to be closer to 0 (than to the 990th percentile) when θ is large and seems to be closer to the 99th percentile for small values of θ.The same behavior (with respect to the second mode) is seen as μ is changing.While it is sort of intuitive to see the behavior of the second mode, it is not easy to explain the logic of having a mode at the origin.
Example 5.2.Here we are interested to see the effect of λ and α 1 on the two measures: the 99th percentile and the mode(s) of the random variable N max .In this example, we fix c = 5, μ = 1.0, θ = 2.0, δ = 1.75, α 2 = 0.3 and calculate the probability mass function of N max .λ is varied in such a way that ρ will take the values 0.2, 0.4, 0.6, and 0.8.Note that as are such that for positive arrivals the average number of customers per batch is a = 2.0 and for the negative arrivals the average is set at d = 3.0.Thus, the parameters for the positive (negative) batch arrivals for these three cases are, respectively, α 1 = 0.5 (α 2 = 2/3), D 1 = 2 (D 2 = 3), and ν 1 = 1.593624 (ν 2 = 2.821439).In Table 5.4, the two measures are displayed for various combinations.
The following observations are obtained for the data seen in Table 5.4.(i) It is interesting to note that the type of distribution considered for the batches of the negative arrivals seems to be some what insensitive to the 99th percentile.However, there seems to be a small dependence when it comes to the distribution of the batch sizes for the regular arrivals.
(ii) With respect to the other measure, there appears to be two modes one at 0 and the other at 6 for all cases.
Example 5.5.The purpose of this example is to see how the probabilities, P(N max = k), behave as functions of ρ.Here we fix c = 5, λ = 3.0, α 1 = 0.5, δ = 1.75, α 2 = 0.3, and θ = 2.0.We vary ρ through changing the values of μ.In Table 5.5, we display some selected probabilities, namely, P(N max = k), k = 0,1,5,10,20.Note that in order to find the number of customers seen in the retrial orbit during a busy period, we do not have to assume the system to be stable.However, as is to be expected, when ρ ≥ 1, the busy period may not terminate (with probability one) and hence the tail of the distributions is heavy.
An examination of the above table indicates the following observations.In the following, let p k = P(N max = k), for k = 0,1,5,10,20.
(i) For both queueing models, the 99th percentile is an increasing function of ρ and a decreasing one of θ.While in general the 99th percentile is lesser for the model with negative arrivals, it is interesting to note that the percentile, for fixed values of ρ, varies in a wider domain for the model with batch arrivals.
(ii) The model with negative arrivals appears to have a second mode when ρ is greater than 0.5.In contrast, the model with batch arrivals exhibits a second mode when the system congestion is high (i.e., when ρ is large and θ is small).
(iii) In both models, for all values of ρ, the percentile seems to approach a limit value as far as the retrial rate increases.The convergence is faster in the model with negative arrivals.