A New Algorithm for the Design of Finite Capacity Service Units

Dario Andrea Bini^{1}^{⋆}, Srinivas Chakravarthy^{2}, and Beatrice Meini^{1}

1 Dipartimento di Matematica Università di Pisa, 56127 Pisa, Italy.

                          bini,meini@dm.unipi.it

2 Department of Science and Mathematics Kettering University, Flint, USA.

                          schakrav@kettering.edu

Abstract. In this paper we present a new algorithm for the design of finite capacity service units with phase type arrivals and phase type services. These systems are very useful in manufacturing and production areas. Using cyclic reduction method that exploits both inner and outer structures of the generator of the Markov process governing PH/PH/1/N queue, we develop a new, fast and numerically efficient algorithm for solving the steady state probability vector. This fast algorithm can interactively be used in solving a number of optimization problems that require searching for an optimal value of N. Comparison of our algorithm to that of block Gauss-Seidel iterative procedure and to the matrix sign function approach for selected examples is presented.

1 Introduction

Finite capacity queueing models play an important role in many areas of engineering, notably in production and manufacturing engineering. Jobs (or messages or customers) arrive to a service facility in a random manner and get processed by the server. In this paper, we consider a single server queueing system in which jobs arrive according to a point process governed by a phase type distribution (PH-distribution) with representation (α, T) of order m. The service times are assumed to be of PH-type with representation (β, S) of order n. The queueing system has a finite capacity of size N. Any arriving job finding the buffer full is considered lost. Thus, we consider a PH/PH/1/N queue. For details on this queueing model, we refer the reader to [5, 14, 15]. Our focus here is to use the latest development in block tridiagonal Toeplitz systems with structured blocks to develop a new, fast and numerically efficient algorithm for solving the steadystate probability vector of the queueing model under study. This algorithm can interactively be used to solve a number of optimization problems useful in practice. A specific optimization problem that is of interest in practice is mentioned in Section 2. The results of the solution of this and other optimization problems will be reported in a paper sequel to this.

⋆ Supported by the MURST project “Advanced numerical methods for scientific computing”.

2 The Model

We consider a single-server queuing system with a finite waiting room of size N. The arrival process is a PH-renewal process where the underlying PH-distribution F (·) has a representation (α, T) of order m. The service times are assumed to be independent, identically distributed with the common PH-distribution G(·) with representation (β, S) of order n.

A probability distribution F (·) on [0, +∞) is a PH-distribution if and only if it is the distribution of the time until absorption in a finite Markov chain of the type

[ \begin{matrix}T & T^{0} \\ 0 & 0 \end{matrix} ]

where T is the transition generator, T^{0} = −T e, e^{T} = (1, 1, …, 1) and 0 denotes the null vector. If α is the initial probability vector, the pair (α, T) is called a representation of F (·). The class of PH-distributions includes many well-known distributions such as generalized Erlang, hyperexponential, etc., as simple special cases and has a number of interesting closure properties. We refer to Neuts [15] for a detailed discussion on this.

Using the PH-formalism as given in Chapter 2 of Neuts [15] efficient algorithms for solving the steady-state probability vectors were developed by Chakravarthy and Neuts [5]. Some interesting relations between various stationary probabilities were derived in that paper. More recently in [9] Latouche and Ramaswami performed a structural analysis in a probabilistic framework that provided a compact representation of the PH/PH/1 queue by means of a QBD process where the block size is reduced to the sum of the number of phases in the arrival and service distributions.

A finite capacity single server queue with phase type arrival and phase type services with a finite buffer size N is denoted by PH/PH/1/N. Any arrivals finding the buffer full will be lost. By keeping track of the phase of the arrivals, the number of customers in the system, and the phase of the service (if the server is busy), this queueing system can be analyzed using matrix formalism [5]. For sake of completeness, we will briefly describe this queueing model. Let W (t), Y (t) and Z(t), respectively, denote the number of customers in the system, the phase of the arrival process, and the phase of the service process (if the server is busy) at time t. Then, {W (t), Y (t), Z(t), t ≥ 0} is a continuous-parameter Markov process with state space given by

Ω = {(0, j), 1 ≤ j ≤ m} ∪ {(i, j, k), 1 ≤ i ≤ N + 1, 1 ≤ j ≤ m, 1 ≤ k ≤ n}.

The Markov process is in state (0, j) if there are no customers in the system and the arrival process is in phase j. The process is in state (i, j, k) if there are i > 0 customers in the system with the arrival process in phase j and the current service is in phase k.

The generator Q of the Markov process is the (m + (N + 1)mn) × (m + (N + 1)mn) matrix of the block-tridiagonal Toeplitz-like form

Q = [ \begin{matrix}T & Ĉ & & & & ◯ \\ B̂ & A & C \\ & B & A & C \\ & & ... & ... & ... \\ & & & B & A & C \\ ◯ & & & & B & K \end{matrix} ]  (1)

where

\begin{matrix}Ĉ = T^{0}α ⊗ β, B̂ = I_{m} ⊗ S^{0}, \\ A = T ⊕ S, C = T^{0}α ⊗ I_{n}, B = I_{m} ⊗ S^{0}β, \\ K = T ⊕ S + T^{0}α ⊗ I_{n}. \end{matrix}  (2)

Here, T^{0} and S^{0} are such that T e+T^{0} = 0 and Se+S^{0} = 0, e = (1, 1, …, 1)^{T}, the symbols ⊗ and ⊕ denote the Kronecker product and the Kronecker sum of matrices, respectively and I_{r} denotes the r × r identity matrix. We refer to Marcus and Minc [10] for details on Kronecker products and sums.

For the sake of clarity it is useful to specify the size of the blocks involved in the matrix Q, specifically, T ∈ IR^{m×m}, S ∈ IR^{n×n}, α ∈ IR^{1×m}, β ∈ IR^{1×n}, T^{0} ∈ IR^{m×1}, S^{0} ∈ IR^{n×1} and the matrix Q in (1) has block size N + 2.

Let y be the steady state probability vector of the state of the Markov process at an arbitrary time. That is, y solves the system yQ = 0 and ye = 1.

Observe that the generator Q has a double structure. More precisely an outer block structure (block tridiagonal block Toeplitz-like), and an inner structure given by the Kronecker sums and products, and by the low rank properties of the blocks B and C.

The main purpose of this paper is to exploit the double structure of the generator in order to develop an ad-hoc algorithm for computing the steady state probability vector y. This task is mandatory if we want to solve the problem for large values of N. In fact, the customary algorithms would impose a huge computational cost. Moreover, the availability of a fast algorithm would improve the efficiency of the interactive solution of the following optimization problem.

Optimization Problem: Find a value of N so that, given real numbers ϵ > 0 and x^{∗} > 0, we have P (X > x^{∗}) < ϵ, where X = X(N) is the stationary waiting time of an admitted customer into the system.

The task of exploiting both the inner and the outer structure of Q is nontrivial, in fact the algorithms (like iterative methods as Gauss-Seidel iteration) that exploit the inner structure do not fully exploit the outer structure, and the algorithms that exploit the outer structure (like cyclic reduction introduced in [4]) may destroy the inner structure.

In practice, for large size of the blocks and for large block size of the matrix Q, it is necessary to exploit both structures. Moreover in solving the optimization problem we need to solve several systems, with different values of the block dimension N + 2.

In the next section we propose a method, based on cyclic reduction, that exploits both inner and outer structures of Q, and that is suitable for solving the above mentioned optimization problem. Specifically, following the results of [13] and [12], due to the block tridiagonal block Toeplitz-like structure of Q, we may express the block components of the invariant vector y in terms of the powers of the matrices R and U, which solve two suitable quadratic matrix equations. The matrices R and U do not depend on the block size N of Q, but only depend on the block entries A, B and C. Thus, in the interactive solution of the optimization problem for different values of N, R and U may be computed only once. In order to compute R and U we apply the cyclic reduction algorithm, where the general iteration formulas are expressed in a different form in order to take advantage of the inner structure.

The computational cost of one step of the new algorithm obtained in this way is reduced from O(m^{3}n^{3}) to O((m + n)^{3}) arithmetic operations.

3 Solving a block tridiagonal Toeplitz system with structured blocks

In this section we are concerned with the design of an efficient numerical method for solving the following linear system

yQ = 0, ye = 1,  (3)

where Q is the matrix of (1).

By partitioning y as y = (y_{i})_{i=0,N+1}, according to the block structure of Q, we may rewrite (3) in terms of a homogeneous linear matrix difference equation

y_{i−1}C + y_{i}A + y_{i+1}B = 0, i = 2, 3, …, N,  (4)

with the boundary conditions

\{ \begin{matrix}y_{0}T + y_{1}B̂ = 0 \\ y_{0}Ĉ + y_{1}A + y_{2}B = 0 \\ y_{N}C + y_{N+1}K = 0. \end{matrix}  (5)

Let the matrices R and U be the minimal nonnegative solutions of

\begin{matrix}C + RA + R^{2}B = 0, \\ U^{2}C + UA + B = 0. \end{matrix}  (6)

Then, by a direct inspection it follows that for any u, v ∈ IR^{1×mn}, the vectors

y_{i} = uR^{i−1} + vU^{N+1−i}, i = 1, …, N + 1,  (7)

solve the equation (4). Furthermore from (5) we have y_{0} = −y_{1}B̂T^{−1}. Imposing the boundary conditions (5) allows us to compute the vectors u, v. In fact, from (7) and (5) we have

[ u v ] [ \begin{matrix}A − B̂T^{−1}Ĉ + RB & R^{N−1}(C + RK) \\ U^{N}(A − B̂T^{−1}Ĉ) + U^{N−1}B & UC + K \end{matrix} ] = 0.  (8)

This system has a solution since its matrix W is singular (it is a simple matter to verify that W e = 0). The same formulas can be found in [1], [12] and [13].

In this way the problem of computing y is reduced to solving the matrix equations (6), the homogeneous system (8) and to applying (7). For the normalization of the vector obtained in this way we may impose the condition

ye = u ∑_{i=0}^{N} R^{i}e + v ∑_{i=0}^{N} U^{i}e − (u + vU^{N})B̂T^{−1}e = 1.

It would be interesting to analyze the sensitivity of the solution y to perturbations of R and U. Even though this problem goes beyond the purpose of this paper, we believe that a conditioning analysis can be carried out by using the tools of [11], more specifically the matrix identity

X^{i} − Y^{i} = ∑_{j=1}^{i} X^{j−1}(X − Y)Y^{i−j}, i ≥ 1,

and its translation in terms of tensor products

vec(X^{i} − Y^{i}) = (∑_{j=1}^{i}(Y^{T})^{i−j} ⊗ X^{j−1}) vec(X − Y),

where vec(A) denotes the vector obtained by arranging the entries of the matrix A columnwise. The conditioning analysis can be completed by applying matrix norms, taking into account that the matrices R and U are stochastic or substochastic.

To solve the matrix equations (6) in [13] the authors apply the logarithmic reduction algorithm of Latouche and Ramaswami [8], in [1] the authors apply the matrix sign function approach of [2]. Here, we apply the algorithm of [4],[3] based on cyclic reduction by fully exploiting the inner structure of A, B and C.

We first recall the following result that reports formula (2.15) of [3].

Theorem 1. The matrices A_{j}, B_{j}, C_{j}, Â_{j}, Ã_{j}, j = 0, 1, …, inductively defined by

\begin{matrix}A_{j+1} = A_{j} − B_{j}A_{j}^{−1}C_{j} − C_{j}A_{j}^{−1}B_{j}, \\ B_{j+1} = −B_{j}A_{j}^{−1}B_{j}, \\ C_{j+1} = −C_{j}A_{j}^{−1}C_{j}, \\ Â_{j+1} = Â_{j} − C_{j}A_{j}^{−1}B_{j}, \\ Ã_{j+1} = Ã_{j} − B_{j}A_{j}^{−1}C_{j}, \end{matrix}  (9)

where A_{0} = Â_{0} = Ã_{0} = A, B_{0} = B, C_{0} = C, are such that det A_{j+1}, det Ã_{j+1}, det Â_{j+1} ≠ 0. Moreover, the matrices R and U solve the matrix equations

\begin{matrix}C + RÂ_{j} + R^{2^{j}+1}B_{j} = 0, \\ U^{2^{j}+1}C_{j} + UÃ_{j} + B = 0, \end{matrix}  (10)

and the matrices R^{2^{j}} and U^{2^{j}} solve the matrix equations

\begin{matrix}C_{j} + R^{2^{j}}A_{j} + R^{2·2^{j}}B_{j} = 0, \\ U^{2·2^{j}}C_{j} + U^{2^{j}}A_{j} + B_{j} = 0. \end{matrix}  (11)

Proof. The matrices U and R satisfy the equations

[ R R^{2} …] [ \begin{matrix}A & C \\ B & A & C \\ & ... & ... & ... \end{matrix} ] = [ −C 0 …] ,

[ U U^{2} …] [ \begin{matrix}A & B \\ C & A & B \\ & ... & ... & ... \end{matrix} ] = [ −B 0 …] ,

respectively. Applying j steps of cyclic reduction to the above systems yields

[ R R^{2^{j}+1} …] [ \begin{matrix}Â_{j} & C_{j} \\ B_{j} & A_{j} & C_{j} \\ & ... & ... & ... \end{matrix} ] = [ −C 0 …] ,

[ U U^{2^{j}+1} …] [ \begin{matrix}Ã_{j} & B_{j} \\ C_{j} & A_{j} & B_{j} \\ & ... & ... & ... \end{matrix} ] = [ −B 0 …] ,

(compare also lemma 3.1 of [4]). By equating the first block component in both sides of the above equations yields (10). Applying the same argument to the systems

[ I R R^{2} …] [ \begin{matrix}A & C \\ B & A & C \\ & ... & ... & ... \end{matrix} ] = [ A + RB 0 …] ,

[ I U U^{2} …] [ \begin{matrix}A & B \\ C & A & B \\ & ... & ... & ... \end{matrix} ] = [ A + UC 0 …] ,

and equating the second block entry at the j-th step of cyclic reduction yield (11). ⊓ ⊔

From the properties of the matrix Q assumed in section 2 and from proposition 1 of [13], it follows that either ρ(U) < 1, or ρ(R) < 1, provided that the matrix UC +A+RB is nonsingular. If ρ(U) < 1 the eigenvalues of U are the nm solutions of the equation p(λ) = det(B + λA + λ^{2}C) = 0 in the open unit disk and the eigenvalues of R are the reciprocal of the remaining solutions (where we assume roots at the infinity if the polynomial p(λ) has degree less than 2nm), hence ρ(R) = 1. If otherwise ρ(R) < 1 the eigenvalues of R are the reciprocal of the nm solutions of the equation p(λ) = 0 of modulus greater than 1 and the eigenvalues of U are the solutions of p(λ) = 0 in the closed unit disk, hence ρ(U) = 1. This fact allows us to prove the following convergence property.

Theorem 2. It holds

\begin{matrix}R = −CÂ_{j}^{−1} + O(r^{2^{j}}), \\ U = −BÃ_{j}^{−1} + O(r^{2^{j}}), \end{matrix}

where r is any number such that min{ρ(U), ρ(R)} < r < 1.

Proof. Without loss of generality assume that ρ(R) < 1. From the first equation of (10) it follows that ||R + CÂ_{j}^{−1}|| = ||R^{2^{j}+1}B_{j}|| for any matrix norm. Now observe that, from the definition of A_{j+1}, since its diagonal entries are nonpositive and B_{j} ≥ 0, C_{j} ≥ 0, it follows that Diag(A_{j+1}) ≥ Diag(A_{j}), where Diag(A) denotes the diagonal matrix having the same diagonal entries of A. This property, together with the relation A_{j}e + B_{j}e + C_{j}e = 0, implies that 0 ≤ B_{j}e + C_{j}e = −A_{j}e ≤ −Diag(A_{j})e ≤ −Diag(A)e. That is B_{j} is bounded by a constant independent of j. From the latter property, since for any ϵ > 0 there exists an operator matrix norm || · || such that ||R|| < ρ(R) + ϵ (compare [7]), it follows the first part of the theorem. For the second part, we observe that ||U^{2^{j}}|| is bounded, in fact ρ(U^{2^{j}}) = 1, so that it is sufficient to prove that ||C_{j}|| = O(r^{2^{j}}) for a suitable operator norm. This follows since from (11) it holds C_{j} = −R^{2^{j}}A_{j} − R^{2·2^{j}}B_{j} and A_{j}, B_{j} have bounded norm. ⊓ ⊔

If UC + A + RB is singular then the conditions ρ(U) < 1 or ρ(R) < 1 may not be satisfied ([13], proposition 1), and the convergence stated in the above theorem is not guaranteed.

From theorem 2 the matrices R and U can be approximated by means of the matrices −CÂ_{j}^{−1}, −BÃ_{j}^{−1}, where j is such that min{||C_{j}||_{∞}, ||B_{j}||_{∞}} < ϵ, for a given input error bound ϵ, where || · ||_{∞} denotes the infinity norm.

In order to compute the blocks A_{j}, B_{j}, C_{j}, Â_{j}, Ã_{j} at each step j, the inner structure of A, B, C can be fully exploited. Indeed, we prove the following result:

Theorem 3. At each step j ≥ 1 of cyclic reduction it holds

\begin{matrix}A_{j} = A − (I_{m} ⊗ S^{0})σ_{j}(α ⊗ I_{n}) − (T^{0} ⊗ I_{n})θ_{j}(I_{m} ⊗ β), \\ B_{j} = −(I_{m} ⊗ S^{0})δ_{j}(I_{m} ⊗ β), \\ C_{j} = −(T^{0} ⊗ I_{n})γ_{j}(α ⊗ I_{n}), \\ Â_{j} = A − (T^{0} ⊗ I_{n})θ_{j}(I_{m} ⊗ β), \\ Ã_{j} = A − (I_{m} ⊗ S^{0})σ_{j}(α ⊗ I_{n}), \end{matrix}  (12)

where δ_{j}, γ_{j}, σ_{j}, θ_{j} are m × m, n × n, m × n, n × m matrices, respectively. The matrices δ_{j}, γ_{j}, σ_{j}, θ_{j}, for j ≥ 2, are defined by the recursions

\begin{matrix}δ_{j+1} = δ_{j} (δ_{1} + [ δ_{1} σ_{1} ] D_{j}V_{j}^{−1} [ \begin{matrix}θ_{1} \\ δ_{1} \end{matrix} ]) δ_{j}, \\ γ_{j+1} = γ_{j} (γ_{1} + [ θ_{1} γ_{1} ] D_{j}V_{j}^{−1} [ \begin{matrix}γ_{1} \\ σ_{1} \end{matrix} ]) γ_{j}, \\ σ_{j+1} = σ_{j} + δ_{j} (σ_{1} + [ δ_{1} σ_{1} ] D_{j}V_{j}^{−1} [ \begin{matrix}γ_{1} \\ σ_{1} \end{matrix} ]) γ_{j}, \\ θ_{j+1} = θ_{j} + γ_{j} (θ_{1} + [ θ_{1} γ_{1} ] D_{j}V_{j}^{−1} [ \begin{matrix}θ_{1} \\ δ_{1} \end{matrix} ]) δ_{j}, \end{matrix}  (13)

for j = 1, 2, …, where

\begin{matrix}δ_{1} = (I_{m} ⊗ β)A^{−1}(I_{m} ⊗ S^{0}), \\ γ_{1} = (α ⊗ I_{n})A^{−1}(T^{0} ⊗ I_{n}), \\ σ_{1} = (I_{m} ⊗ β)A^{−1}(T^{0} ⊗ I_{n}), \\ θ_{1} = (α ⊗ I_{n})A^{−1}(I_{m} ⊗ S^{0}), \end{matrix}  (14)

and

D_{j} = [ \begin{matrix}σ_{j} & 0 \\ 0 & θ_{j} \end{matrix} ] , V_{j} = I_{m+n} − [ \begin{matrix}θ_{1} & γ_{1} \\ δ_{1} & σ_{1} \end{matrix} ] D_{j}.

Proof. Let us observe that, for the properties of Kronecker product we have B = I_{m} ⊗ S^{0}β = (I_{m} ⊗ S^{0})(I_{m} ⊗ β) and C = T^{0}α ⊗ I_{n} = (T^{0} ⊗ I_{n})(α ⊗ I_{n}). Due to this structure, by applying one step of cyclic reduction by means of (9) with A_{0} = A, B_{0} = B, C_{0} = C we arrive at the formulas (14) for representing A_{1}, B_{1}, C_{1}, Â_{1} and Ã_{1} in terms of δ_{1}, γ_{1}, σ_{1} and θ_{1}. For the remaining part, we proceed by induction on j. Let us assume that (12) hold for j and prove them for j +1. We need to describe the structure of A_{j}^{−1}. Let us write A_{j} = A−Z_{j}W_{j}, where

Z_{j} = [ I_{m} ⊗ S^{0} T^{0} ⊗ I_{n} ] D_{j}, W_{j} = [ \frac{α ⊗ I_{n}}{I_{m} ⊗ β} ] ,

and observe that V_{j} = I − W_{j}A^{−1}Z_{j}. From the nonsingularity of A and A_{j} (compare theorem 1), since det(A_{j}) = det(A) det(V_{j}), it follows that det(V_{j}) ≠ 0. Whence, applying the Sherman Woodbury Morrison formula A_{j}^{−1} = A^{−1} + A^{−1}Z_{j}(I − W_{j}A^{−1}Z_{j})^{−1}W_{j}A^{−1}, we obtain that

A_{j}^{−1} = A^{−1} + A^{−1} [ I_{m} ⊗ S^{0} T^{0} ⊗ I_{n} ] D_{j}V_{j}^{−1} [ \frac{α ⊗ I_{n}}{I_{m} ⊗ β} ] A^{−1}.  (15)

By applying the equations (9), in the light of the inductive hypothesis, we may complete the proof. We show this for the matrix B_{j+1}, we leave the details on the remaining part to the reader. We have

B_{j+1} = −B_{j}A_{j}^{−1}B_{j} = −(I_{m} ⊗ S^{0}){δ_{j}(I_{m} ⊗ β)A_{j}^{−1}(I_{m} ⊗ S^{0})δ_{j}}(I_{m} ⊗ β).

Whence B_{j+1} has the structure (12) where δ_{j+1} = δ_{j}(I_{m} ⊗ β)A_{j}^{−1}(I_{m} ⊗ S^{0})δ_{j}. By replacing A_{j}^{−1} with (15) and by using the definitions of (14), we obtain the first equation of (13). ⊓ ⊔

In this way, the computation of the mn × mn blocks A_{j}, B_{j}, C_{j}, Ã_{j}, Â_{j}, is reduced to computing the blocks of lower dimensions δ_{j}, γ_{j}, σ_{j}, θ_{j}. Thus the computational cost of one step is reduced from O(m^{3}n^{3}) to O((m + n)^{3}) arithmetic operations.

It is interesting to point out that the matrices δ_{1}, γ_{1}, σ_{1}, θ_{1} of (14) are the same matrices displayed in formulas (33)–(36) of [9], obtained by means of probabilistic arguments.

Also the matrices R and U are structured. Indeed, it is easy to verify that

\begin{matrix}Â_{j}^{−1} = A^{−1} + A^{−1}(T^{0} ⊗ I_{n})θ_{j}(I_{m} − σ_{1}θ_{j})^{−1}(I_{m} ⊗ β)A^{−1}, \\ Ã_{j}^{−1} = A^{−1} + A^{−1}(I_{m} ⊗ S^{0})σ_{j}(I_{n} − θ_{1}σ_{j})^{−1}(α ⊗ I_{n})A^{−1}. \end{matrix}

Hence, R and U are approximated by the matrices R_{j} = −CÂ_{j}^{−1}, U_{j} = −BÃ_{j}^{−1}, respectively, for an enough large value of j, such that

\begin{matrix}R_{j} = −(T^{0} ⊗ I_{n}) ((α ⊗ I_{n}) + γ_{1}θ_{j}(I_{m} − σ_{1}θ_{j})^{−1}(I_{m} ⊗ β)) A^{−1}, \\ U_{j} = −(I_{m} ⊗ S^{0}) ((I_{m} ⊗ β) + δ_{1}σ_{j}(I_{n} − θ_{1}σ_{j})^{−1}(α ⊗ I_{n})) A^{−1}. \end{matrix}  (16)

Thus the matrices R and U have rank m and n, respectively, and

\begin{matrix}R_{j}^{k} = −(T^{0} ⊗ I_{n})ν_{R}^{k−1} ((α ⊗ I_{n}) + γ_{1}θ_{j}(I_{m} − σ_{1}θ_{j})^{−1}(I_{m} ⊗ β)) A^{−1}, \\ U_{j}^{k} = −(I_{m} ⊗ S^{0})ν_{U}^{k−1} ((I_{m} ⊗ β) + δ_{1}σ_{j}(I_{n} − θ_{1}σ_{j})^{−1}(α ⊗ I_{n})) A^{−1}, k ≥ 1, \end{matrix}

where ν_{R} and ν_{U} are n × n and m × m matrices, respectively, given by

\begin{matrix}ν_{R} = − (γ_{1} + γ_{1}θ_{j}(I_{m} − σ_{1}θ_{j})^{−1}σ_{1}) , \\ ν_{U} = − (δ_{1} + δ_{1}σ_{j}(I_{n} − θ_{1}σ_{j})^{−1}θ_{1}) . \end{matrix}  (17)

From (16) it follows that the matrices R and U have a tensor structure as stated by the following

Theorem 4. The solutions R and U of (6) are such that

\begin{matrix}R = −(T^{0} ⊗ I_{n}) ((α ⊗ I_{n}) + γ_{1}θ^{∗}(I_{m} − σ_{1}θ^{∗})^{−1}(I_{m} ⊗ β)) A^{−1} \\ U = −(I_{m} ⊗ S^{0}) ((I_{m} ⊗ β) + δ_{1}σ^{∗}(I_{n} − θ_{1}σ^{∗})^{−1}(α ⊗ I_{n})) A^{−1}. \end{matrix}

where θ^{∗} = lim_{j} θ_{j}, σ^{∗} = lim_{j} σ_{j} and θ_{j}, σ_{j} are defined in (13).

Proof. The limits lim_{j} θ_{j}, lim_{j} σ_{j} exist since the sequences θ_{j} and σ_{j} are monotone and bounded. The expressions for R and U follow directly from (16). ⊓ ⊔

The structure of R and U can thus be exploited in the computation of the entries of the matrix in system (8), and in the computation of y_{i}, by means of formula (7).

The resulting algorithm is resumed below:

Algorithm 1. Computing R and U.

Input: Positive integers n, m, the vectors α ∈ IR^{m}, β ∈ IR^{n}, the matrices T ∈ IR^{m×m}, S ∈ IR^{n×n}, and an error bound ϵ > 0.

Output: Approximations θ̃ ∈ IR^{n×m}, σ̃ ∈ IR^{m×n} to θ^{∗}, σ^{∗} providing approximations R̃, Ũ to the solutions R and U of the equation (6) by means of theorem 4.

Computation: 1. Compute T^{0} = −T e, S^{0} = −Se.

2. Compute A = T ⊕ S and A^{−1}.

3. Compute the matrices δ_{1}, γ_{1}, σ_{1}, θ_{1} by means of (14).

4. For j = 1, 2, … compute the matrices δ_{j}, γ_{j}, σ_{j}, θ_{j} by means of (13) until

min{||γ_{j}||_{∞}, ||δ_{j}||_{∞}} < ϵ.

5. Set θ̃ = θ_{j}, σ̃ = σ_{j}.

For the computation of the vector y we rely on the following:

Algorithm 2. Computing the vector y.

Input: Positive integer n, m, N, the vectors α ∈ IR^{m}, β ∈ IR^{n}, the matrices

T ∈ IR^{m×m}, S ∈ IR^{n×n}, θ^{∗} ∈ IR^{n×m}, σ^{∗} ∈ IR^{m×n}, the matrices δ_{1}, γ_{1}, σ_{1}, θ_{1}. Output: Approximations of the components y_{0}, …, y_{N+1} of the steady state

vector. Computation: 1. Apply the repeated squaring technique for computing ν_{R}^{N−2}

and ν_{U}^{N−2} in ⌈log_{2}(N − 2)⌉ steps (compare (17)).

2. Solve the system (8).

3. By means of successive vector-matrix products, compute uR^{i} and vU^{i} for i = 0, …, N and compute y_{i} by means of (7) for i = 1, …, N + 1. For this computation exploit the tensor product structure of U and R given by theorem 4.

4. Compute y_{0} = −y_{1}B̂T^{−1}.

5. Normalize the vector y.

The computational costs of algorithms 1 and 2 can easily be estimated. The cost of algorithm 1, dominated by the computation of a matrix inversion of size m + n and a finite number of matrix products of size m + n, amounts to O((m+n)^{3} log log ϵ^{−1}) arithmetic operations (ops) due to the double exponential convergence of the scheme. Concerning algorithm 2, the repeated squaring stage costs O((max{m, n})^{3} log N) ops, the computation of (8) and its solution costs O((m+n)^{3}) ops, the computation of y_{i} by means of (7) costs O((max{m, n})^{2}N) ops. Thus, the overall computation of algorithm 2 costs O((m+n)^{3} log N +(m+ n)^{2}N) ops.

The steady state vector could be directly computed by applying the folding algorithm of [16] or the cyclic reduction method as described in [3]. However this approach would require a number of steps in algorithm 1 proportional to log N and would involve a large memory storage needed for the back substitution stage. In our approach the number of steps, proportional to log log ϵ^{−1}, does not depend on N, moreover we do not have to store the matrices obtained in the intermediate steps of the computation.

For PH/PH/1 queues, where the generator matrix is infinite and has the same structure of (1), the components of the steady state vector are given by y_{i} = wR^{i}, i = 1, 2, …, where w is a suitable row vector and R is the solution of the first matrix equation in (6). Thus the approach proposed in this paper can be used for computing the matrix R and the steady state vector y.

4 Numerical Results

We implemented our algorithm using Fortran 90 and compared it with block Gauss-Seidel iterative algorithm, for the computation of the vector y. We also performed a comparison with the matrix sign function algorithm for the computation of the matrices R and U by using the software TELPACK that can

be taken at the url http://www.cstp.umkc.edu/org/tn/telpack/home.html.

The latter software does not exploit the inner structure of the problem; the specific structure of the matrix Q is exploited when using the Gauss-Seidel method.

The software implementing our algorithm can be downloaded at the url

http://www.dm.unipi.it/pages/meini/public html.

We ran our benchmarks on an alpha platform with the 53 bit standard IEEE arithmetic.

We considered the following test problems, where the arrivals and services are taken to be Erlang or mixed Erlang or hyperexponential. As is known, the convergence of the Gauss-Seidel procedure is very slow, especially for the hyperexponential case, and hence it is imperative to have a procedure with faster convergence in all cases.

Problem 1. n = m = 5, α = β = (1, 0, 0, 0, 0), T = H(10/3), S = H(5), where

H(r) = [ \begin{matrix}−r & r & 0 & 0 & 0 \\ 0 & −r & r & 0 & 0 \\ 0 & 0 & −r & r & 0 \\ 0 & 0 & 0 & −r & r \\ 0 & 0 & 0 & 0 & −r \end{matrix} ] .

Problem 2. n = m = 5, α = (1, 0, 0, 0, 0), β = (0.7, 0, 0, 0.3, 0), T = H(10/3),

S = [ \begin{matrix}−8 & 8 & 0 & 0 & 0 \\ 0 & −8 & 8 & 0 & 0 \\ 0 & 0 & −8 & 0 & 0 \\ 0 & 0 & 0 & −0.8135593 & 0.8135593 \\ 0 & 0 & 0 & 0 & −0.8135593 \end{matrix} ] .

Problem 3. m = 5, n = 2, β = (1, 0, 0, 0, 0), α = (0.9, 0.1), S = H(5),

T = [ \begin{matrix}−2 & 0 \\ 0 & −0.095238 \end{matrix} ] .

Problem 4. m = 5, n = 2, α = (0.9, 0.1), β = (0.7, 0, 0, 0.3, 0),

\begin{matrix}S = [ \begin{matrix}−8 & 8 & 0 & 0 & 0 \\ 0 & −8 & 8 & 0 & 0 \\ 0 & 0 & −8 & 0 & 0 \\ 0 & 0 & 0 & −0.8135593 & 0.8135593 \\ 0 & 0 & 0 & 0 & −0.8135593 \end{matrix} ] , \\ T = [ \begin{matrix}−2 & 0 \\ 0 & −0.095238 \end{matrix} ] . \end{matrix}

CR  GS
N Time Iter. Res. Time Iter Res.
50 0.02 5 6.6e-16 0.71 287 2.9e-15
100 0.02 5 6.6e-16 1.85 430 3.5e-15
200 0.26 5 6.6e-16 5.84 693 3.3e-15
400 0.48 5 6.6e-16 20.19 1181 4.3e-15
800 0.85 5 6.6e-16 73.5 2115 3.7e-15
1600 1.61 5 6.6e-16 270.80 3891 4.1e-15

Table 1. Problem 1

CR  GS
N Time Iter. Res. Time Iter Res.
50 0.01 7 3.0e-16 1.45 662 2.4e-15
100 0.01 7 3.0e-16 3.62 850 2.3e-15
200 0.02 7 3.0e-16 9.89 1173 2.6e-15
400 0.41 7 3.0e-16 29.81 1752 2.5e-15
800 0.65 7 3.0e-16 96.8 2794 2.4e-15
1600 1.26 7 3.0e-16 328.76 4755 2.5e-15

Table 2. Problem 2

The CPU time in seconds, the number of iterations and the residual errors of the Gauss-Seidel (GS) algorithm and our algorithm (CR) in the computation of y for several values of N, are displayed in Tables 1 through 4.

Examining the entries in these four tables, it is evident that our algorithm performs far better than Gauss-Seidel algorithm in terms of CPU time and the number of iterations. In fact, the speed-up factor is greater than 20 for all the values of N and increases as N grows. For problem 4, with the value N = 1600, the speed-up factor is about 382. We observe that, despite the residual errors have the same order of magnitude, the numerical behavior of the two algorithms is different. For example, if we look at the components of the vector y provided by the algorithms, we see that the new algorithm generates an absolute error, that is about the machine precision, uniformly distributed in all the components of y. This means that components less than the machine precision may be affected by large relative errors. Whereas the Gauss-Seidel iteration seems to provide components affected by a relative error roughly bounded by the machine precision.

If the user is interested in an approximation with a uniform relative error, it is more convenient to apply a suitable number of steps of the Gauss-Seidel iteration to the vector provided by our algorithm in order to refine the smallest components. This kind of refinement does not affect much the computational time of the whole procedure and provides a very accurate solution.

Table 5 reports the CPU time in seconds and the residual error of the computation of the matrices R and U performed by means of the cyclic reduction (CR) and the matrix sign function (MSF) approach implemented in the package TELPACK of [2].

CR  GS
N Time Iter. Res. Time Iter Res.
50 0.05 9 2.0e-16 1.02 1275 1.3e-16
100 0.06 9 2.0e-16 3.75 2416 1.5e-16
200 0.10 9 2.0e-16 10.17 3318 1.2e-16
400 0.15 9 2.0e-16 25.75 4201 1.3e-16
800 0.27 9 2.0e-16 73.86 5591 1.4e-16
1600 0.52 9 2.0e-16 199.02 8060 1.4e-16

Table 3. Problem 3

CR  GS
N Time Iter. Res. Time Iter Res.
50 0.06 9 3.2e-17 0.94 1180 1.6e-16
100 0.06 9 3.2e-17 3.77 2449 1.6e-16
200 0.10 9 3.2e-17 10.85 3551 1.8e-16
400 0.16 9 3.2e-17 27.81 4525 1.6e-16
800 0.30 9 3.2e-17 73.89 5976 1.6e-16
1600 0.59 9 3.2e-17 210.14 8491 1.7e-16

Table 4. Problem 4

CR  MSF
Time Res. Time Res.
Problem 1 0.06 1.1e-16 0.38 1.4e-14
Problem 2 0.05 1.1e-16 0.42 8.8e-14
Problem 3 0.04 2.2e-16 0.12 1.5e-14
Problem 4 0.04 1.1e-16 0.13 1.6e-14

Table 5. CPU time and residual for the computation of R and U.

Besides being faster than the matrix sign function approach, our method provides a better approximation with a smaller residual error.

Acknowledgments. We wish to thank the anonymous referees for their comments and suggestions that improved the presentation.

References

1. N. Akar and K. Sohraby. Finite and infinite QBD chains: a simple unifying algorithmic approach. In Proc. IEEE Infocom 97. 1997.

2. N. Akar and K. Sohraby. An invariant subspace approach in M/G/1 and G/M/1 type Markov chains. Commun. Statist. Stochastic Models, 13:381–416, 1997.

3. D. A. Bini and B. Meini. On cyclic reduction applied to a class of Toeplitz-like matrices arising in queueing problems. In Computations with Markov Chains, Ed. W. J. Stewart, Kluwer Academic Publisher, 21–38, 1995.

4. D. A. Bini and B. Meini. On the solution of a nonlinear matrix equation arising in queueing problems. SIAM J. Matrix Anal. Appl., 17:906–926, 1996.

5. S. Chakravarthy and M. F. Neuts. Algorithms for the Design of Finite-Capacity Service Units, Naval Research Logistics, 36, 147-165, 1989.

6. S. D. Conte and Carl de Boor. Numerical Analysis: An Algorithmic Approach, McGraw-Hill, New York, 1972.

7. G. H. Golub and C. F. Van Loan. Matrix Computations. The Johns Hopkins University Press, Baltimore, 1989.

8. G. Latouche and V. Ramaswami. A logarithmic reduction algorithm for Quasi-Birth-Death processes. J. Appl. Probability, 30:650–674, 1993.

9. G. Latouche and V. Ramaswami. The PH/PH/1 queue at epochs of queue size change. Queueing Syst. Appl., 25:97–114, 1997.

10. M. Marcus and H. Minc. A Survey of Matrix Theory and Matrix Inequalities, Allyn and Bacon, Boston, MA, 1964.

11. B. Meini. New convergence results on functional iteration techniques for the numerical solution of M/G/1 type Markov chains. Numer. Math., 78:39–58, 1997.

12. V. Naoumov. Matrix-multiplicative approach to quasi-birth-and-death processes analysis. In Matrix-analytic methods in stochastic models, Eds. S.R. Chakravarthy and A.S. Alfa, Marcel Dekker, Inc., New York, 87–106, 1996.

13. V. Naoumov, U. R. Krieger, and D. Wagner. Analysis of a multi-server delay-loss system with a general markovian arrival process. In Matrix-analytic methods in stochastic models, Eds. S.R. Chakravarthy and A.S. Alfa, Marcel Dekker, Inc., New York, 43–66, 1996.

14. M. F. Neuts. Matrix-Geometric Solutions in Stochastic Models: An algorithmic Approach, The Johns Hopkins University Press, Baltimore, MD, 1981.

15. M. F. Neuts. Structured Stochastic Matrices of M/G/1 Type and Their Applications, Marcel Dekker Inc., New York, 1989.

16. J. Ye and San-Qi Li. Analysis of multi-media traffic queues with finite buffer and overload control–Part I: algorithm. In Proc. IEEE Infocom 91, Bal Harbour, FL, pages 1464–1474. 1991.