A sequential update algorithm for computing the stationary distribution vector in upper block-Hessenberg Markov chains

TL;DR

This paper introduces a new sequential update algorithm for accurately computing the stationary distribution in upper block-Hessenberg Markov chains, leveraging LBCL-augmented truncations and linear fractional programming.

Contribution

It develops a novel algorithm using LBCL-augmented truncations and LFP problems to compute the exact stationary distribution without approximation.

Findings

Algorithm effectively computes stationary distributions for complex queues.

Demonstrated applicability to BMAP/M/∞ and M/M/s retrial queues.

Provides bounds and convergence guarantees for the method.

Abstract

This paper proposes a new algorithm for computing the stationary distribution vector in continuous-time upper block-Hessenberg Markov chains. To this end, we consider the last-block-column-linearly-augmented (LBCL-augmented) truncation of the (infinitesimal) generator of the upper block-Hessenberg Markov chain. The LBCL-augmented truncation is a linearly-augmented truncation such that the augmentation distribution has its probability mass only on the last block column. We first derive an upper bound for the total variation distance between the respective stationary distribution vectors of the original generator and its LBCL-augmented truncation. Based on the upper bound, we then establish a series of linear fractional programming (LFP) problems to obtain augmentation distribution vectors such that the bound converges to zero. Using the optimal solutions of the LFP problems, we construct a matrix-infinite-product (MIP) form of the original (i.e., not approximate) stationary distribution vector and develop a sequential update algorithm for computing the MIP form. Finally, we demonstrate the applicability of our algorithm to BMAP/M/$\infty$ queues and M/M/$s$ retrial queues.

Equations

$$\mathbb{S}$$

$${\mathbb{L}}_k$$

$$\bm{Q}=\bordermatrix{&\mathbb{L}_{0}&\mathbb{L}_{1}&\mathbb{L}_{2}&\mathbb{L}_{3}&\cdots\cr\mathbb{L}_{0}&\bm{Q}_{0,0}&\bm{Q}_{0,1}&\bm{Q}_{0,2}&\bm{Q}_{0,3}&\cdots\cr\mathbb{L}_{1}&\bm{Q}_{1,0}&\bm{Q}_{1,1}&\bm{Q}_{1,2}&\bm{Q}_{1,3}&\cdots\cr\mathbb{L}_{2}&\bm{O}&\bm{Q}_{2,1}&\bm{Q}_{2,2}&\bm{Q}_{2,3}&\cdots\cr\mathbb{L}_{3}&\bm{O}&\bm{O}&\bm{Q}_{3,2}&\bm{Q}_{3,3}&\cdots\cr~{}\vdots&\vdots&\vdots&\vdots&\vdots&\ddots}.$$

$$\mathbf{\pi }$$

$$\displaystyle\hskip 0.50003pt{}_{(n)}\bm{Q}=\left(\begin{array}[]{ccccccc}\bm{Q}_{0,0}&\bm{Q}_{0,1}&\bm{Q}_{0,2}&\cdots&\bm{Q}_{0,n-2}&\bm{Q}_{0,n-1}&\bm{Q}_{0,n}\\ \bm{Q}_{1,0}&\bm{Q}_{1,1}&\bm{Q}_{1,2}&\cdots&\bm{Q}_{1,n-2}&\bm{Q}_{1,n-1}&\bm{Q}_{1,n}\\ \bm{O}&\bm{Q}_{2,1}&\bm{Q}_{2,2}&\cdots&\bm{Q}_{2,n-2}&\bm{Q}_{2,n-1}&\bm{Q}_{2,n}\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\ \bm{O}&\bm{O}&\bm{O}&\cdots&\bm{Q}_{n-1,n-2}&\bm{Q}_{n-1,n-1}&\bm{Q}_{n-1,n}\\ \bm{O}&\bm{O}&\bm{O}&\cdots&\bm{O}&\bm{Q}_{n,n-1}&\bm{Q}_{n,n}\\ \end{array}\right).$$

$${}_{(n)}\overline{\mathbf{Q}}={}_{(n)}\mathbf{Q}-{}_{(n)}\mathbf{Q}\mathbf{e}{}_{(n)}\mathbf{\alpha },n\in {\mathbb{Z}}_{+},$$

$${}_{(n)}\overline{\mathbf{\pi }}=\frac{{}_{(n)}\mathbf{\alpha }(-{}_{(n)}\mathbf{Q}{)}^{-1}}{{}_{(n)}\mathbf{\alpha }(-{}_{(n)}\mathbf{Q}{)}^{-1}\mathbf{e}},n\in {\mathbb{Z}}_{+},$$

$${}_{(n)}\overline{\mathbf{\pi }}{}_{(n)}\overline{\mathbf{Q}}=0,{}_{(n)}\overline{\mathbf{\pi }}\ge 0,{}_{(n)}\overline{\mathbf{\pi }}\mathbf{e}=1;$$

$$\hskip 0.50003pt{}_{(n)}\widehat{\bm{\alpha}}=\bordermatrix{&\mathbb{S}_{n-1}&\mathbb{L}_{n}\cr&\bm{0}&\bm{\alpha}_{n}}.$$

$${}_{(n)}\hat{\mathbf{Q}}={}_{(n)}\mathbf{Q}-{}_{(n)}\mathbf{Q}\mathbf{e}{}_{(n)}\hat{\mathbf{\alpha }},n\in {\mathbb{Z}}_{+}.$$

$${}_{(n)}\hat{\mathbf{\pi }}=\frac{{}_{(n)}\hat{\mathbf{\alpha }}(-{}_{(n)}\mathbf{Q}{)}^{-1}}{{}_{(n)}\hat{\mathbf{\alpha }}(-{}_{(n)}\mathbf{Q}{)}^{-1}\mathbf{e}},n\in {\mathbb{Z}}_{+}.$$

$${\mathbf{\pi }}^{(N)}=\frac{({\mathbf{\pi }}_0,{\mathbf{\pi }}_1,\dots ,{\mathbf{\pi }}_N)}{{\sum }_{\ell =0}^N{\mathbf{\pi }}_{\ell }},$$

$${}_{(n)}\hat{\mathbf{\pi }}|{}_{(n)}\hat{\mathbf{Q}}-\mathbf{Q}|$$

$$\hskip 0.50003pt{}_{(n)}\widehat{\bm{\alpha}}\bm{v}=\bordermatrix{&\mathbb{S}_{n}&\overline{\mathbb{S}}_{n}\cr&\hskip 0.50003pt{}_{(n)}\widehat{\bm{\alpha}}&\bm{0}}\bm{v}=\bordermatrix{&\mathbb{S}_{n-1}&\mathbb{L}_{n}&\overline{\mathbb{S}}_{n}\cr&\bm{0}&\bm{\alpha}_{n}&\bm{0}}\bm{v}=\bm{\alpha}_{n}\bm{v}_{n},$$

$${}_{(n)}\hat{\mathbf{\pi }}$$

$$(-{}_{(n)}\mathbf{Q}{)}^{-1}$$

$${}_{(n)}\hat{\mathbf{\alpha }}(-{}_{(n)}\mathbf{Q}{)}^{-1}={\mathbf{\alpha }}_n({}_{(n)}{\mathbf{X}}_{n,0},{}_{(n)}{\mathbf{X}}_{n,1},\dots ,{}_{(n)}{\mathbf{X}}_{n,n}),n\in {\mathbb{Z}}_{+}.$$

$${}_{(n)}\hat{\mathbf{\pi }}=\frac{{\mathbf{\alpha }}_n({}_{(n)}{\mathbf{X}}_{n,0},{}_{(n)}{\mathbf{X}}_{n,1},\dots ,{}_{(n)}{\mathbf{X}}_{n,n})}{{\mathbf{\alpha }}_n{\sum }_{\ell =0}^n{}_{(n)}{\mathbf{X}}_{n,\ell }\mathbf{e}},n\in {\mathbb{Z}}_{+},$$

$${}_{(n)}{\hat{\mathbf{\pi }}}_k=\frac{{\mathbf{\alpha }}_n{}_{(n)}{\mathbf{X}}_{n,k}}{{\mathbf{\alpha }}_n{\sum }_{\ell =0}^n{}_{(n)}{\mathbf{X}}_{n,\ell }\mathbf{e}},n\in {\mathbb{Z}}_{+},k\in {\mathbb{Z}}_n.$$

$$\sum _{\ell =0}^n{}_{(n)}{\mathbf{X}}_{n,\ell }\mathbf{e}$$

$$\hskip 0.50003pt{}_{(n)}\bm{Q}=\left(\begin{array}[]{cccc|c}&&&&\,\bm{Q}_{0,n}\\ &&\hskip 0.50003pt{}_{(n-1)}\bm{Q}&&\,\vdots\\ \rule[-5.69054pt]{0.0pt}{11.38109pt}&&&&\,\bm{Q}_{n-1,n}\\ \hline\cr\rule{0.0pt}{11.38109pt}\bm{O}&\cdots&\bm{O}&\bm{Q}_{n,n-1}&\,\bm{Q}_{n,n}\end{array}\right),\qquad n\in\mathbb{N}.$$

$${}_{(n)}{\mathbf{X}}_{n,n}$$

$$=$$

$${}_{(n)}{\mathbf{X}}_{n,k}$$

$${\mathbf{U}}_n^{\ast }$$

$${\mathbf{U}}_{n,k}$$

$$\bm{U}_{n}^{*}=\left\{\begin{array}[]{l@{~~~}l}(-\bm{Q}_{0,0})^{-1},\hfil~{}~{}~{}&n=0,\\ \left(-\bm{Q}_{n,n}-\displaystyle\sum_{\ell=0}^{n-1}\bm{U}_{n,\ell}\bm{Q}_{\ell,n}\right)^{-1},\hfil~{}~{}~{}&n\in\mathbb{N}.\end{array}\right.$$

$${}_{(n)}{\mathbf{X}}_{n,k}={\mathbf{U}}_n^{\ast }{\mathbf{U}}_{n,k},k\in {\mathbb{Z}}_n,$$

$$\bm{U}_{n,k}=\left\{\begin{array}[]{l@{~~}l}(\bm{Q}_{n,n-1}\bm{U}_{n-1}^{*})(\bm{Q}_{n-1,n-2}\bm{U}_{n-2}^{*})\cdots(\bm{Q}_{k+1,k}\bm{U}_{k}^{*}),\hfil~{}~{}&k\in\mathbb{Z}_{n-1},\\ \bm{I},\hfil~{}~{}&k=n.\end{array}\right.$$

$${\mathbf{U}}_{n,k}={\mathbf{Q}}_{n,n-1}{\mathbf{U}}_{n-1}^{\ast }{\mathbf{U}}_{n-1,k},k\in {\mathbb{Z}}_{n-1},$$

$${\mathbf{\pi }}_{\ell }={\mathbf{\pi }}_n{\mathbf{U}}_{n,\ell },n\in \mathbb{N},\ell \in {\mathbb{Z}}_n.$$

$${}_{(n)}{\hat{\mathbf{\pi }}}_k$$

Full text

A sequential update algorithm for computing

the stationary distribution vector in upper block-Hessenberg Markov chains111Published online in Queueing Systems on February 21, 2019 (doi: 10.1007/s11134-019-09599-x)

---

Hiroyuki Masuyama222E-mail: masuyama@sys.i.kyoto-u.ac.jp

Department of Systems Science, Graduate School of Informatics, Kyoto University

Kyoto 606-8501, Japan

Abstract

[TABLE]

[TABLE]

1 Introduction

This paper considers an upper block-Hessenberg Markov chain in continuous time. To describe such a Markov chain, we first introduce some symbols. Let $\mathbb{R}_{+}$ denote the set of all nonnegative real numbers, i.e., $\mathbb{R}_{+}=[0,\infty)$. Let $\mathbb{N}=\{1,2,3,\dots\}$, $\mathbb{Z}_{+}=\{0,1,2,\dots\}$, and $\mathbb{Z}_{n}=\{0,1,\dots,n\}$ for $n\in\mathbb{Z}_{+}$. We then introduce some sets of pairs of integers:

[TABLE]

where $\mathbb{M}_{k}=\{1,2,\dots,M_{k}\}\subset\mathbb{N}$. We also define $(k,i;\ell,j)$ as an ordered pair $((k,i),(\ell,j))$ in $\mathbb{S}^{2}$. Furthermore, we define $\bm{e}=(1,1,\dots)^{\top}$, which has an appropriate (finite or infinite) number of ones.

Let $\{(X(t),J(t));t\in\mathbb{R}_{+}\}$ denote a regular-jump bivariate Markov chain with state space $\mathbb{S}$ (see [3, Chapter 8, Definition 2.5] for the definition of regular-jump Markov chains). Let $\bm{Q}:=(q(k,i;\ell,j))_{(k,i;\ell,j)\in\mathbb{S}^{2}}$ denote the (infinitesimal) generator of the Markov chain $\{(X(t),J(t))\}$, which is in an upper block-Hessenberg form:

[TABLE]

We refer to $\{(X(t),J(t))\}$ as the upper block-Hessenberg Markov chain (which may be called the level-dependent M/G/1-type Markov chain) and refer to $X(t)$ and $J(t)$ as the level variable and the phase variable, respectively. Note that if $(X(t),J(t))\in\mathbb{L}_{k}$ then $X(t)=k$ and thus $\mathbb{L}_{k}$ is called level $k$.

Throughout the paper, unless otherwise stated, we assume that $\{(X(t),J(t))\}$ is ergodic (i.e., irreducible, aperiodic and positive recurrent). We then define $\bm{\pi}:=(\pi(k,i))_{(k,i)\in\mathbb{S}}>\bm{0}$ as the unique stationary distribution vector of the ergodic generator $\bm{Q}$ (see, e.g., [1, Chapter 5, Theorems 4.4 and 4.5]). By definition, $\bm{\pi}\bm{Q}=\bm{0}$ and $\bm{\pi}\bm{e}=1$. For later use, we also define $\bm{\pi}_{k}=(\pi(k,i))_{i\in\mathbb{M}_{k}}$ for $k\in\mathbb{Z}_{+}$ and partition $\bm{\pi}$ as

[TABLE]

It is, in general, difficult to obtain an explicit expression of $\bm{\pi}=(\bm{\pi}_{0},\bm{\pi}_{1},\dots)$. Thus, we study the computation of the stationary distribution vector $\bm{\pi}$ through a linearly augmented truncation of the ergodic generator $\bm{Q}$. The linearly augmented truncation is described below.

Let $\hskip 0.50003pt{}_{(n)}\bm{Q}:=(\hskip 0.50003pt{}_{(n)}q(k,i;\ell,j))_{(k,i;\ell,j)\in(\mathbb{S}_{n})^{2}}$, $n\in\mathbb{Z}_{+}$, denote the northwest corner truncation of the ergodic generator $\bm{Q}$, which is given by

[TABLE]

We then define $\hskip 0.50003pt{}_{(n)}\overline{\bm{Q}}:=(\hskip 0.50003pt{}_{(n)}\overline{q}(k,i;\ell,j))_{(k,i;\ell,j)\in(\mathbb{S}_{n})^{2}}$, $n\in\mathbb{Z}_{+}$, as a $Q$-matrix (diagonally dominant matrix with nonpositive diagonal elements and nonnegative off-diagonal ones; see, e.g., [1, Section 2.1]) such that

[TABLE]

where $\hskip 0.50003pt{}_{(n)}\bm{\alpha}:=(\hskip 0.50003pt{}_{(n)}\alpha(k,i))_{(k,i)\in\mathbb{S}_{n}}$ is a probability vector. We refer to $\hskip 0.50003pt{}_{(n)}\overline{\bm{Q}}$ as the linearly augmented truncation of $\bm{Q}$. We also refer to $\hskip 0.50003pt{}_{(n)}\bm{\alpha}$ as the augmentation distribution vector.

Let $\hskip 0.50003pt{}_{(n)}\overline{\bm{\pi}}=(\hskip 0.50003pt{}_{(n)}\overline{\pi}(k,i))_{(k,i)\in\mathbb{S}_{n}}$, $n\in\mathbb{Z}_{+}$, denote

[TABLE]

where $(-\hskip 0.50003pt{}_{(n)}\bm{Q})^{-1}$ exists due to the ergodicity of $\bm{Q}$. From (1.3) and (1.4), we have

[TABLE]

that is, $\hskip 0.50003pt{}_{(n)}\overline{\bm{\pi}}$ is a stationary distribution vector of the linearly augmented truncation $\hskip 0.50003pt{}_{(n)}\overline{\bm{Q}}$. Furthermore, as $n\to\infty$, each element of $\hskip 0.50003pt{}_{(n)}\overline{\bm{Q}}$ converges to the corresponding one of $\bm{Q}$. Thus, we can expect $\hskip 0.50003pt{}_{(n)}\overline{\bm{\pi}}$ to be an approximation to $\bm{\pi}$. This is why we refer to $\hskip 0.50003pt{}_{(n)}\overline{\bm{\pi}}$ as the linearly augmented truncation approximation to $\bm{\pi}$.

We note that if the augmentation distribution vector $\hskip 0.50003pt{}_{(n)}\bm{\alpha}$ has its probability mass only on the last block (i.e., $\mathbb{L}_{n}$) then the linearly augmented truncation $\hskip 0.50003pt{}_{(n)}\overline{\bm{Q}}$ inherits upper block-Hessenberg structure from the original generator $\bm{Q}$. To utilize this tractable structure, we focus on a special linearly augmented truncation $\hskip 0.50003pt{}_{(n)}\overline{\bm{Q}}$ with $\hskip 0.50003pt{}_{(n)}\bm{\alpha}=\hskip 0.50003pt{}_{(n)}\widehat{\bm{\alpha}}$, where $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\alpha}}$ is a probability vector such that

[TABLE]

For convenience, we refer to such a linearly augmented truncation as a last-block-column-linearly-augmented truncation (LBCL-augmented truncation).

We now define $\hskip 0.50003pt{}_{(n)}\widehat{\bm{Q}}:=(\hskip 0.50003pt{}_{(n)}\widehat{q}(k,i;\ell,j)_{(k,i;\ell,j)\in(\mathbb{S}_{n})^{2}}$, $n\in\mathbb{Z}_{+}$, as the LBCL-augmented truncation of $\bm{Q}$, that is, a $Q$-matrix such that

[TABLE]

We also define $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}:=(\hskip 0.50003pt{}_{(n)}\widehat{\pi}(k,i))_{(k,i)\in\mathbb{S}_{n}}$, $n\in\mathbb{Z}_{+}$, as

[TABLE]

Note that $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$ is equal to $\hskip 0.50003pt{}_{(n)}\overline{\bm{\pi}}$ in (1.4) with $\hskip 0.50003pt{}_{(n)}\bm{\alpha}=\hskip 0.50003pt{}_{(n)}\widehat{\bm{\alpha}}$; that is, $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$ is a stationary distribution vector of the LBCL-augmented truncation $\hskip 0.50003pt{}_{(n)}\widehat{\bm{Q}}$. Hence, we call $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$ the last-block-column-linearly-augmented truncation approximation (LBCL-augmented truncation approximation) to $\bm{\pi}$.

In this paper, we propose a new algorithm for computing the original stationary distribution vector $\bm{\pi}$ by using the LBCL-augmented truncation approximation $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$. In fact, $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$ does not necessarily converge to $\bm{\pi}$ as $n\to\infty$ (see Section 2.3). We solve such a problem by choosing $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\alpha}}$ adaptively for each $n\in\mathbb{Z}_{+}$. To achieve this, we first derive an upper bound for the total variation distance between $\bm{\pi}$ and $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$. With this upper bound, we establish a series of linear fractional programming (LFP) problems for finding $\{\hskip 0.50003pt{}_{(n)}\widehat{\bm{\alpha}};n\in\mathbb{Z}_{+}\}$ such that $\{\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}};n\in\mathbb{Z}_{+}\}$ converges to $\bm{\pi}$. Fortunately, the optimal solutions of the LFP problems are explicitly obtained. Thus, we can readily construct a convergent sequence of LBCL-augmented truncation approximations, which yields a matrix-infinite-product (MIP) form of $\bm{\pi}$. We note that the LFP problems are not given in advance but are formulated successively while constructing the MIP form. As a result, we can develop a sequential update algorithm for computing $\bm{\pi}$.

We now review related work. Some researchers [2, 4, 21] have studied the computation of level-dependent quasi-birth-and-death processes (LD-QBDs), which belong to a special case of upper block-Hessenberg Markov chains. These previous studies propose algorithms for computing the conditional stationary distribution vector $\bm{\pi}^{(N)}$:

[TABLE]

where $N\in\mathbb{N}$ is the truncation parameter that should be determined so that $\bm{\pi}^{(N)}$ is sufficiently close to $\bm{\pi}$. Takine [24] develops an algorithm for computing $\bm{\pi}^{(N)}$ of a special upper block-Hessenberg Markov chain, which assumes that, for all sufficiently large $n\in\mathbb{Z}_{+}$, the $\bm{Q}_{n,n-1}$ are nonsingular and the $\bm{Q}_{n,n}$ are of the same order (see Assumption 1 therein). These additional assumptions in [24] are removed by Kimura and Takine [10]. Besides, Shin and Pearce [23], Li et al. [16], and Klimenok and Dudin [11] modify transition rates (or transition probabilities) such that they are eventually level independent, and then these researchers establish algorithms for computing approximately the stationary distribution vectors of upper block-Hessenberg Markov chains.

The algorithms proposed in [4, 23] have update procedures to improve their outputs, like our algorithm. However, their update procedures need to recompute, from scratch, most components of their new outputs every time. On the other hand, our algorithm utilizes the components of the current result, together with some additional computation, to generate an updated result. This is a remarkable feature of our algorithm.

The rest of this paper is divided into four sections. Section 2 describes preliminary results on the LBCL-augmented truncation approximation for upper block-Hessenberg Markov chains. Section 3 proposes a sequential update algorithm that generates a sequence of LBCL-augmented truncation approximations converging to the original stationary distribution vector. Section 4 demonstrates the applicability of the proposed algorithm. Finally, Section 5 provides concluding remarks.

2 The LBCL-augmented truncation approximation

This section consists of three subsections. In Section 2.1, we show a matrix-product form of the LBCL-augmented truncation approximation $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$. In Section 2.2, we derive an error bound for $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$, more specifically, an upper bound for the total variation distance between $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$ and $\bm{\pi}$. In Section 2.3, we provide an example such that $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$ does not converge to $\bm{\pi}$ as $n\to\infty$.

Before entering the body of this section, we describe our notation. For any matrix $\bm{M}$ (resp. vector $\bm{m}$), let $|\bm{M}|$ (resp. $|\bm{m}|$) denote the matrix (resp. vector) obtained by taking the absolute value of each element of $\bm{M}$ (resp. $\bm{m}$). A finite matrix is treated, if necessary, as an infinite matrix that keeps the existing elements in their original positions and has an infinite number of zeros in the other positions. Such treatment is also applied to finite vectors. Thus, for example, it follows from (1.6) that

[TABLE]

where $\hskip 0.50003pt{}_{(n)}\bm{Q}_{>n}=(q(k,i;\ell,j))_{(k,i;\ell,j)\times\mathbb{S}_{n}\times\overline{\mathbb{S}}_{n}}$. It also follows from (1.5) that, for any column vector $\bm{v}:=(v(k,i))_{(k,i)\in\mathbb{S}}$,

[TABLE]

where $\bm{v}_{n}=(v(k,i))_{(k,i)\in\mathbb{L}_{n}}=(v(n,i))_{i\in\mathbb{M}_{n}}$ for $n\in\mathbb{Z}_{+}$. Furthermore, we use the following notation: If a sequence $\{\bm{Z}_{n};n\in\mathbb{Z}_{+}\}$ of finite matrices (or vectors) converges element-wise to an infinite matrix (or vector) $\bm{Z}$, then we denote this convergence by $\lim_{n\to\infty}\bm{Z}_{n}=\bm{Z}$. We also define the empty sum as zero (e.g., $\sum_{k=1}^{0}\cdot=0$).

2.1 A matrix-product form

We partition $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$ and $(-\hskip 0.50003pt{}_{(n)}\bm{Q})^{-1}$ level-wise as follows:

[TABLE]

From (2.6) and (1.5), we have

[TABLE]

Substituting (2.7) into (1.7) yields

[TABLE]

which leads to

[TABLE]

Note that, because $(-\hskip 0.50003pt{}_{(n)}\bm{Q})^{-1}\hskip 0.50003pt{}_{(n)}\bm{Q}=-\bm{I}$, the inverse matrix $(-\hskip 0.50003pt{}_{(n)}\bm{Q})^{-1}\geq\bm{O}$ has no zero rows. Therefore,

[TABLE]

We derive a matrix-product form of $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}_{k}$, $k\in\mathbb{Z}_{n}$, from (2.8). To do this, we need some preparation. We first partition $\hskip 0.50003pt{}_{(n)}\bm{Q}$ as

[TABLE]

From this equation and (2.6), we have the following (see the last two equations in [7, Section 0.7.3]): For $n\in\mathbb{N}$,

[TABLE]

and

[TABLE]

We also define

[TABLE]

for $n\in\mathbb{Z}_{+}$. It then follows from (2.15), (2.17), and (2.20) that

[TABLE]

Using $\bm{U}_{n}^{*}$ and $\bm{U}_{n,k}$, we can express $\hskip 0.50003pt{}_{(n)}\bm{X}_{n,k}$, $k\in\mathbb{Z}_{n}$, as follows.

Lemma 2.1

For $n\in\mathbb{Z}_{+}$,

[TABLE]

and

[TABLE]

Proof.

Combining (2.16) with (2.17) and (2.20), we have (2.22). Furthermore, applying (2.22) to (2.20) yields

[TABLE]

which leads to (2.23). ∎∎

Remark 2.1

A result similar to Lemma 2.1 is presented in Shin [22] under the condition that $\hskip 0.50003pt{}_{(n)}\bm{Q}$ is block tridiagonal (see Theorem 2.1 therein).

Remark 2.2

The matrices $\hskip 0.50003pt{}_{(n)}\bm{X}_{k,\ell}$, $\ell\in\mathbb{Z}_{k}$, and $\bm{U}_{n,\ell}$, $\ell\in\mathbb{Z}_{n-1}$, have probabilistic interpretations. The $(i,j)$-th element of $\hskip 0.50003pt{}_{(n)}\bm{X}_{k,\ell}$ represents the expected total sojourn time in state $(\ell,j)$ before the first visit to $\overline{\mathbb{S}}_{n}$ (i.e., to any state above level $n$) starting from state $(k,i)$ (see, e.g., [13, Theorem 2.4.3]). Furthermore, the $(i,j)$-th element of $\bm{U}_{n,\ell}$ represents the expected total sojourn time in state $(\ell,j)$ before the first visit to $\overline{\mathbb{S}}_{n}$ starting from state $(n,i)$, measured per unit of time spent in state $(n,i)$. Thus, we have (see [13, Equation (5.33)])

[TABLE]

We now obtain a matrix-product form of $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}_{k}$, $k\in\mathbb{Z}_{n}$, by substituting (2.22) into (2.8).

Lemma 2.2

[TABLE]

Remark 2.3

Equations (2.9) and (2.22) lead to

[TABLE]

2.2 An error bound

In this subsection, we present an error bound for the LBCL-augmented truncation approximation $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$ to $\bm{\pi}$. The error bound is used to develop an algorithm for computing $\bm{\pi}$ in the next section.

To derive the error bound, we assume a Foster-Lyapunov drift condition.

Condition 1

The generator $\bm{Q}$ is irreducible, and there exist a constant $b\in(0,\infty)$, a finite set $\mathbb{C}\subset\mathbb{S}$, and a positive column vector $\bm{v}:=(v(k,i))_{(k,i)\in\mathbb{S}}$ such that $\inf_{(k,i)\in\mathbb{S}}v(k,i)>0$ and

[TABLE]

where $\bm{1}_{\mathbb{B}}:=(1_{\mathbb{B}}(k,i))_{(k,i)\in\mathbb{S}}$, $\mathbb{B}\subseteq\mathbb{S}$, denotes a column vector defined by

[TABLE]

Remark 2.4

Recall that $\bm{Q}$ is the generator of the regular-jump Markov chain $\{(X(t),J(t))\}$ (see Section 1) and thus $\bm{Q}$ is stable, i.e., $|q(\ell,j;\ell,j)|<\infty$ for all $(\ell,j)\in\mathbb{S}$ (see, e.g., [3, Chapter 8, Definition 2.4 and Theorem 3.4]). The irreducibility of $\bm{Q}$ and the finiteness of $\mathbb{C}$ imply that

[TABLE]

which shows that $\mathbb{C}$ is a small set (see, e.g., [12]). Therefore, if Condition 1 holds, then the irreducible generator $\bm{Q}$ is ergodic (see, e.g., [12, Theorem 1.1]).

Let $\bm{\Phi}^{(\beta)}:=(\phi^{(\beta)}(k,i;\ell,j)_{(k,i;\ell,j)\in\mathbb{S}}$ denote a stochastic matrix such that

[TABLE]

where $\bm{P}^{(t)}:=(p^{(t)}(k,i;\ell,j))_{(k,i;\ell,j)\in\mathbb{S}}$, $t\in\mathbb{R}_{+}$, is the transition matrix function of the Markov chain $\{(X(t),J(t))\}$ with generator $\bm{Q}$, i.e.,

[TABLE]

Because $\bm{Q}$ is ergodic, we have $\bm{\Phi}^{(\beta)}>\bm{O}$. We also define $\overline{\phi}_{\mathbb{C}}^{(\beta)}$, $\beta>0$, as

[TABLE]

We then have the following result from [17, Theorem 2.1] with $\bm{f}=\bm{g}=\bm{e}$.

Proposition 2.1

Under Condition 1, the following holds:

[TABLE]

where, for any vector $\bm{m}:=(m(i))$, $\|\bm{m}\|$ denotes the total variation norm of $\bm{m}$, i.e., $\|\bm{m}\|=\sum_{i}|m(i)|$.

From Proposition 2.1, we derive a more informative bound for $\|\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}-\bm{\pi}\|$. For this purpose, we define $\bm{U}^{*}_{n,k}$, $n\in\mathbb{Z}_{+}$, $k\in\mathbb{Z}_{n}$, as

[TABLE]

We also define $\bm{u}^{*}_{n}:=(u_{n}^{*}(i))_{i\in\mathbb{M}_{n}}$, $n\in\mathbb{Z}_{+}$, as

[TABLE]

where $\bm{u}^{*}_{n}>\bm{0}$ due to (2.26). Using (2.29) and (2.30), we rewrite (2.25) as

[TABLE]

Theorem 2.1

If Condition 1 holds, then

[TABLE]

where $E(\,\cdot\,):=E^{(\beta)}(\,\cdot\,)$, called the error bound function, is given by

[TABLE]

with $\beta>0$.

Remark 2.5

The error bound function $E$ has a free parameter $\beta>0$ involved in the intractable factor $\overline{\phi}_{\mathbb{C}}^{(\beta)}$. Thus, it is, in general, difficult to discuss theoretically how $\beta$ impacts on the decay speed of $E$. Through numerical experiments, Masuyama [20] investigates such a problem for the last-column block-augmented truncation, though the function $E$ is referred to therein as the error decay function, instead of the error bound function. Note that the last-column block-augmented truncation belongs to the class of block-augmented truncations (see [15] for details). Therefore, the last-column block-augmented truncation is indeed different from our LBCL-augmented truncation, though they are fairly similar.

*Proof of Theorem 2.1 * Suppose that Condition 1 holds. It then follows from (2.4) that

[TABLE]

Substituting (2.39) into (2.28), we have, for $n\in\mathbb{Z}_{+}$,

[TABLE]

where the last equality holds due to (2.5) and $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\alpha}}\bm{e}=1$ for $n\in\mathbb{Z}_{+}$. Because $\bm{Q}\bm{e}=\bm{0}$,

[TABLE]

Incorporating this into (2.44), we have

[TABLE]

Note here that (1.7) and $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\alpha}}\bm{e}=1$ yield

[TABLE]

Thus, we can rewrite (2.45) as

[TABLE]

Furthermore, from (2.7), (2.22), and (2.30), we have

[TABLE]

Applying this equation and (2.31) to (2.46), we obtain

[TABLE]

which results in (2.32) together with (2.33). ∎

2.3 A counterexample to convergence

In the previous subsection, we have established the error bound for the LBCL-augmented truncation approximation $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$. We note that, even if the truncation parameter $n$ goes to infinity, $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$ does not necessarily converge to $\bm{\pi}$, in general. However, it always holds that $\lim_{n\to\infty}\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}=\bm{\pi}$ for special upper block-Hessenberg Markov chains such that the block matrices $\bm{Q}_{k,\ell}$ are scalars. For such a special case, Gibson and Seneta [6] prove that any augmented truncation approximation converges to the original stationary distribution as the truncation parameter goes to infinity (see Theorem 2.2 therein). Of course, this is not the case for general upper block-Hessenberg Markov chains. Indeed, we introduce a counterexample [9].

Fix $\mathbb{M}_{n}=\{1,2\}$ for all $n\in\mathbb{Z}_{+}$, and assume that the block matrices $\bm{Q}_{k,\ell}$ satisfy the following:

[TABLE]

where the symbol “  $\star$  ” denotes some nonzero element. In this case, $\bm{Q}$ is irreducible (see Figure 1), but $\mathbb{S}_{2k-1}$ is not reachable from state $(2k,2)$ avoiding $\overline{\mathbb{S}}_{2k}$.

Thus, the probabilistic interpretations (see Remark 2.2) of the matrices $\bm{U}_{2k}^{*}=\hskip 0.50003pt{}_{(2k)}\bm{X}_{2k,2k}$ and $\bm{U}_{2k,\ell}$, $\ell\in\mathbb{Z}_{2k-1}$, implies that

[TABLE]

We now assume that $\bm{Q}$ is ergodic. We then set

[TABLE]

which implies that $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}$ is the last-column-augmented truncation approximation to the stationary distribution vector $\bm{\pi}>\bm{0}$ of $\bm{Q}$. Applying (2.59), (2.62), and (2.63) to (2.25) yields

[TABLE]

which shows that $\{\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}};n\in\mathbb{Z}_{+}\}$ does not converge to $\bm{\pi}$ in the present setting.

The example presented here implies that, in some cases, the convergence of $\{\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}\}$ to $\bm{\pi}$ can require an adaptive choice of the augmentation distribution vector $\hskip 0.50003pt{}_{(n)}\bm{\alpha}$, depending on $n$. We discuss this problem in the next section.

3 Main results

This section is divided into three subsections. In Section 3.1, we formulate linear fractional programming (LFP) problems for finding augmentation distribution vectors such that the error bound function $E$ converges to zero, i.e., $\lim_{n\to\infty}E(n)=0$. In Section 3.2, using the optimal solutions of these LFP problems, we construct an MIP form of $\bm{\pi}$. In Section 3.3, we present a sequential update algorithm for computing the MIP form.

In this section, we assume that Condition 1 holds, as in Section 2.2. We also assume that $n$ takes an arbitrary value in $\mathbb{Z}_{+}$, unless otherwise stated.

3.1 LFP problems for an MIP form of the stationary distribution vector

Consider the following LFP problem for each $n\in\mathbb{Z}_{+}$:

[TABLE]

It follows from (2.33), (3.1a), and (3.1d) that

[TABLE]

Furthermore, let $\bm{\alpha}_{n}^{*}:=(\alpha_{n}^{*}(j))_{j\in\mathbb{M}_{n}}$ denote a probability vector such that

[TABLE]

where

[TABLE]

We then have the following theorem.

Theorem 3.1

For each $n\in\mathbb{Z}_{+}$, the probability vector $\bm{\alpha}_{n}^{*}$ is an optimal solution of the LFP problem (3.1).

Proof.

From (3.3) and (3.4), we have

[TABLE]

which leads to $\bm{y}_{n}\geq\xi_{n}\bm{u}_{n}^{*}>\bm{0}$. Thus, for any $1\times M_{n}$ probability vector $\bm{p}_{n}$, we obtain

[TABLE]

Therefore, $\bm{\alpha}_{n}^{*}$ is an optimal solution of the LFP problem (3.1). ∎∎

3.2 An MIP form of the stationary distribution vector

Let $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}^{*}:=(\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}_{0}^{*},\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}_{1}^{*},\dots,\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}_{n}^{*})$ denote a probability vector such that

[TABLE]

where ${\rm row}\{\,\cdot\,\}_{j}$ denotes the $j$-th row of the matrix in the brackets. Note here that $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}_{k}^{*}$ is equal to $\hskip 0.50003pt{}_{(n)}\widehat{\bm{\pi}}_{k}$ in (2.31) with $\bm{\alpha}_{n}=\bm{\alpha}_{n}^{*}$. Therefore, it follows from Theorem 2.1 that

[TABLE]

where function $E^{*}$ is equal to $E$ given in (3.2) with $\bm{\alpha}_{n}=\bm{\alpha}_{n}^{*}$; that is,

[TABLE]

To proceed further, we assume the following.

Condition 2

[TABLE]

where $\bm{\Delta}_{n}:=(\Delta_{n}(i,j))_{i,j\in\mathbb{M}_{n}}$ denotes an $M_{n}\times M_{n}$ diagonal matrix such that

[TABLE]

Lemma 3.1

Suppose that Conditions 1 and 2 hold. We then have

[TABLE]

Remark 3.1

If $\bm{Q}$ is bounded, i.e., $\sup_{(n,i)\in\mathbb{S}}\Delta_{n}(i,i)<\infty$, then Condition 2 is reduced to

[TABLE]

*Proof of Lemma 3.1 * To prove this lemma, we require the following proposition (which is proved in Appendix A).

Proposition 3.1

Under Condition 1,

[TABLE]

Let $\widetilde{\bm{\alpha}}_{n}$ denote

[TABLE]

which is well-defined due to Proposition 3.1. Note that $\widetilde{\bm{\alpha}}_{n}$ is a feasible solution of the LFP problem (3.1). Thus, by the optimality of $\bm{\alpha}_{n}^{*}$, we have

[TABLE]

It follows from (3.1a) and (3.12) that

[TABLE]

where the second equality holds due to (3.1d). It also follows from (2.24), (2.29), and (2.30) that

[TABLE]

Substituting these equations into (3.13), and using (3.11), we obtain

[TABLE]

Consequently, the proof of (3.10) is completed by showing that the right-hand side of (3.14) converges to zero as $n\to\infty$.

It follows from (2.27) that, for all $n\in\mathbb{Z}_{+}$ and $k\in\mathbb{Z}_{n}$,

[TABLE]

and thus

[TABLE]

where the last inequality holds due to (3.8). Therefore, by the dominated convergence theorem,

[TABLE]

It also follows from (3.8) that

[TABLE]

Combining (3.15), (3.16), and $\sum_{\ell=0}^{\infty}\bm{\pi}_{\ell}\bm{e}=1$, we obtain

[TABLE]

which completes the proof. ∎

The following theorem is a consequence of Lemma 3.1 together with (3.6) and (3.7).

Theorem 3.2

Suppose that Conditions 1 and 2 hold. We then have

[TABLE]

and thus (3.6) yields

[TABLE]

Proof.

We prove only (3.17). It follows from (3.1d), (3.3), and $\inf_{(k,i)\in\mathbb{S}}v(k,i)>0$ (see Condition 1) that

[TABLE]

Therefore, (3.10) implies that

[TABLE]

which yields

[TABLE]

Applying (3.19) and Lemma 3.1 to (3.7) results in (3.17). ∎∎

Theorem 3.2 yields a matrix-infinite-product (MIP) form of $\bm{\pi}=(\bm{\pi}_{0},\bm{\pi}_{1},\dots)$ under Conditions 1 and 2. This is summarized in the following corollary.

Corollary 3.1

If Conditions 1 and 2 hold, then

[TABLE]

or equivalently,

[TABLE]

where

[TABLE]

Proof.

Suppose that Conditions 1 and 2 hold. It then follows from (3.5) and (3.18) that

[TABLE]

which shows that (3.20) holds. Furthermore, combining (2.29) with (2.23) and (3.22) yields, for $n\in\mathbb{Z}_{+}$,

[TABLE]

Using this and (2.30), we can rewrite (3.20) as (3.21).∎∎

Remark 3.2

Theorem 3.2 ensures that the convergence in (3.20) and (3.21) is uniform for $k\in\mathbb{Z}_{+}$.

Remark 3.3

Another MIP form of $\bm{\pi}_{k}$ is presented in the preprint [19], under some technical conditions different from Conditions 1 and 2.

3.3 A sequential update algorithm for the MIP form

In this subsection, we propose an algorithm for computing $\bm{\pi}$, based on Theorem 3.2 and Corollary 3.1. Our algorithm sequentially updates the LBCL-augmented truncation approximation so that it converges to the MIP form (3.20) of $\bm{\pi}$.

To efficiently achieve this update procedure, we derive recursive formulas. Combining (2.29) with (2.21) and (2.23), we have

[TABLE]

Using (2.30), (3.23a), and (3.23d), we also obtain

[TABLE]

Furthermore, (2.23) and (2.29) yield

[TABLE]

Substituting this into (2.21) leads to

[TABLE]

Our algorithm is composed of the equations (3.23)–(3.25), Theorem 3.2, and Corollary 3.1.

Remark 3.4

Equation (3.18) leads to

[TABLE]

Therefore, our algorithm iterates Step 4 only a finite number of times.

Remark 3.5

Step (4.b) computes $\bm{U}_{n}^{*}$ by (3.25). The $(i,j)$-th element of $\bm{U}_{n}^{*}$ is the expected total sojourn time in state $(n,j)$ before the first visit to $\overline{\mathbb{S}}_{n}$ starting from state $(n,i)$. Thus, $\bm{T}_{n}^{*}=(-\bm{U}_{n}^{*})^{-1}$, defined in (A.1), is a non-conservative $Q$-matrix that governs the transient transitions of an absorbing Markov chain obtained by observing $\{(X(t),J(t))\}$ when it is in $\mathbb{L}_{n}$ during the first passage time to $\overline{\mathbb{S}}_{n}$ starting from $\mathbb{L}_{n}$. This consideration indicates $\bm{U}_{n}^{*}=(-\bm{T}_{n}^{*})^{-1}$ can be efficiently computed (see [14, Proposition 1]), provided that $\bm{T}_{n}^{*}$ is given.

Remark 3.6

Generally, our algorithm computes the infinite sum $\sum_{\ell=n+1}^{\infty}\bm{Q}_{k,\ell}\bm{v}_{\ell}$ to obtain $\bm{y}_{n}$ in (3.1d). However, this infinite sum can be calculated in many practical cases associated with queueing models (as implied by the examples in the next section). Moreover, if $\bm{Q}$ is an LD-QBD generator, or equivalently, $\bm{Q}_{k,\ell}=\bm{O}$ for $k\in\mathbb{Z}_{+}$ and $|\ell-k|>1$, then $\bm{y}_{n}$ is expressed without any infinite sum:

[TABLE]

Furthermore, a noteworthy fact is that computing the infinite sum $\sum_{\ell=n+1}^{\infty}\bm{Q}_{k,\ell}\bm{v}_{\ell}$ is not always necessary even if $\bm{Q}$ is not an LD-QBD generator. To demonstrate this, suppose that we have an explicit expression for $\bm{w}_{k,n}$, $k,n\in\mathbb{Z}_{+}$, such that

[TABLE]

It then follows from (3.14) and (3.16) that

[TABLE]

Thus, we modify Step (4.d.i) as follows: Compute

[TABLE]

and find

[TABLE]

Despite this modification, our update algorithm works well.

4 Applicability of the proposed algorithm

This section demonstrates the applicability of our algorithm. To this end, we consider a BMAP/M/$\infty$ queue and M/M/$s$ retrial queue, respectively, in Sections 4.1 and 4.2. For each model, we present a sufficient condition for Conditions 1 and 2, under which our update algorithm works well.

4.1 BMAP/M/$\infty$ queue

This subsection considers a BMAP/M/$\infty$ queue. The system has an infinite number of servers. Customers arrive at the system according to a batch Markovian arrival process (BMAP) (see, e.g., [18]). Arriving customers are immediately served, and their service times are independent and identically distributed (i.i.d.) with an exponential distribution having mean $\mu^{-1}$.

Let $\{N(t);t\in\mathbb{R}_{+}\}$ denote the counting process of arrivals from the BMAP; that is, $N(t)$ is equal to the total number of arrivals during the time interval $[0,t]$, where $N(0)=0$. Let $\{J(t);t\in\mathbb{R}_{+}\}$ denote the background Markov chain of the BMAP, which is defined on state space $\mathbb{M}=\{1,2,\dots,M\}\subset\mathbb{N}$. We assume that the bivariate stochastic process $\{(N(t),J(t));t\in\mathbb{R}_{+}\}$ is a continuous-time Markov chain which follows the transition law given by

[TABLE]

where $a(t)=o(b(t))$ represents $\lim_{t\to 0}a(t)/b(t)=0$. Thus, the BMAP is characterized by $\{\bm{D}_{n};n\in\mathbb{Z}_{+}\}$, where $\bm{D}_{n}=(D_{n,i,j})_{i,j\in\mathbb{M}}$ for $n\in\mathbb{Z}_{+}$. Moreover, $\bm{D}:=\sum_{n\in\mathbb{Z}_{+}}\bm{D}_{n}$ is the generator of the background Markov chain $\{J(t);t\in\mathbb{R}_{+}\}$. As usual, we assume that $\bm{D}$ is irreducible and

[TABLE]

Let $X(t)$, $t\in\mathbb{R}_{+}$, denote the number of customers in the system at time $t$. It then follows that $\{(X(t),J(t));t\in\mathbb{R}_{+}\}$ is a continuous-time Markov chain on state space $\mathbb{S}:=\mathbb{Z}_{+}\times\mathbb{M}$ with generator $\bm{Q}$ given by

[TABLE]

where $\mathbb{L}_{k}=\{k\}\times\mathbb{M}$ (i.e., $M_{k}=M$) for all $k\in\mathbb{Z}_{+}$ and

[TABLE]

We now suppose that, for some $C>0$,

[TABLE]

and let

[TABLE]

where “${\rm e}$” denotes Napier’s constant. Clearly, $\sum_{k=1}^{\infty}\log(k+{\rm e})\bm{D}_{k}\bm{e}<C\bm{e}$ due to (4.7). Thus, Condition 1 holds for generator $\bm{Q}$ in (4.2) (see [25, Lemma 1]).

It remains to verify that Condition 2 holds. From (3.9) and (4.2), we have $\Delta_{k}(i,i)=|q(k,i;k,i)|=k\mu+|D_{0,i,i}|$ and thus Condition 2 is reduced to

[TABLE]

Therefore, we show that (4.8) holds.

We begin with the following lemma.

Lemma 4.1

Let $V$ denote a function on $\mathbb{R}_{+}$ such that

[TABLE]

If (4.7) holds, then there exist some $K\in\mathbb{Z}_{+}$ and $\theta>0$ such that

[TABLE]

Proof.

Because $\lim_{x\to\infty}V(x)=\infty$, it suffices to prove that

[TABLE]

It follows from (4.6) that, for $k\in\mathbb{N}$,

[TABLE]

Furthermore, $V$ is differentiable and convex. Thus, we have

[TABLE]

Using this inequality and (4.9), we obtain

[TABLE]

Applying (4.13) to (4.12) yields

[TABLE]

and therefore (4.11) holds if

[TABLE]

Consequently, our goal is to prove (4.14).

We note that $V\geq 1$ is log-concave, which implies the following: For any $x,y\in\mathbb{R}_{+}$ such that $x+y>0$,

[TABLE]

These inequalities yield

[TABLE]

which leads to

[TABLE]

Using (4.15) and (4.9), we obtain, for all $k\in\mathbb{Z}_{+}$,

[TABLE]

where the last inequality is due to (4.7). Thus, by the dominated convergence theorem and (4.9), we obtain

[TABLE]

which shows that (4.14) holds. ∎∎

Let $\widetilde{\bm{v}}:=(\widetilde{v}(k,i))_{(k,i)\in\mathbb{Z}_{+}\times\mathbb{M}}$ and $\widetilde{\bm{f}}:=(\widetilde{f}(k,i))_{(k,i)\in\mathbb{Z}_{+}\times\mathbb{M}}$ denote column vectors such that

[TABLE]

where $K\in\mathbb{Z}_{+}$ and $\theta>0$ satisfying (4.10). It then follows from Lemma 4.1 that, for some $\widetilde{b}>0$,

[TABLE]

which yields $\bm{\pi}\widetilde{\bm{f}}<\widetilde{b}$. Combining this inequality and (4.18) results in (4.8). We have confirmed that Condition 2 is satisfied. As a result, our algorithm is always applicable to BMAP/M/$\infty$ queues satisfying (4.7).

4.2 M/M/$s$ retrial queue

In this subsection, we consider an M/M/$s$ retrial queue (which is sometimes called an M/M/$s$/$s$ retrial queue). The system has $s$ ($s\in\mathbb{N}$) servers but no real waiting room. Primary customers (which originate from the exterior) arrive to the system according to a Poisson process with rate $\lambda\in(0,\infty)$. If an arriving primary customer finds an idle server, then the customer occupies the server, otherwise it joins the orbit (i.e., the virtual waiting room). Customers in the orbit are referred to as retrial customers. Each retrial customer stays in the orbit for an exponentially distributed time with mean $\eta^{-1}\in(0,\infty)$, independently of all the other events. After the sojourn in the orbit, a retrial customer tries to occupy one of idle servers. If such a retrial customer finds no idle servers, then it goes back to the orbit; that is, becomes a retrial customer again. We assume that the service times of primary and retrial customers are i.i.d. with an exponential distribution having mean $\mu^{-1}\in(0,\infty)$.

Let $X(t)$, $t\in\mathbb{R}_{+}$, denote the number of customers in the orbit at time $t$. Let $J(t)$, $t\in\mathbb{R}_{+}$, denote the number of busy servers at time $t$. The stochastic process $\{(X(t),J(t));t\in\mathbb{R}_{+}\}$ is a level-dependent quasi-birth-and-death process (LD-QBD) on state space $\mathbb{S}:=\mathbb{Z}_{+}\times\mathbb{Z}_{s}$ with generator $\bm{Q}$ given by

[TABLE]

where $\mathbb{L}_{k}=\{k\}\times\mathbb{Z}_{s}$ for $k\in\mathbb{Z}_{+}$, and where

[TABLE]

and

[TABLE]

with

[TABLE]

We now assume that the stability condition $\rho:=\lambda/(s\mu)<1$ holds. It then follows that the LD-QBD $\{(X(t),J(t))\}$ is ergodic (see, e.g., [5, Section 2.2]) and thus has the unique stationary distribution vector $\bm{\pi}=(\pi(k,i))_{(k,i)\in\mathbb{S}}$. Under this stability condition, we show that Conditions 1 and 2 are satisfied, which requires the following proposition.

Proposition 4.1** **([20, Lemma 4.1])

Suppose that $\rho=\lambda/(s\mu)<1$. For $k\in\mathbb{Z}_{+}$, let $\bm{v}_{k}=(v(k,j))_{j\in\mathbb{Z}_{s}}$ be given by

[TABLE]

where $\alpha$, $\gamma$, and $c$ are positive constants such that

[TABLE]

Furthermore, let

[TABLE]

where $x\vee y=\max(x,y)$. Under these conditions, the generator $\bm{Q}$ of the LD-QBD, characterized by (4.19)–(4.38), satisfies

[TABLE]

We note that $c\bm{v}\geq\bm{e}.$ Proposition 4.1 thus shows that Condition 1 is satisfied. Moreover, Theorem 1 in [8] states that, for a certain constant $c_{0}>0$,

[TABLE]

where $a_{1}(x)\stackrel{{\scriptstyle x}}{{\sim}}a_{2}(x)$ represents $\lim_{x\to\infty}a_{1}(x)/a_{2}(x)=1$. Combining (4.41)–(4.43) yields

[TABLE]

which implies that Condition 2 is satisfied. Consequently, our algorithm is always applicable to stable M/M/$s$ retrial queues.

5 Concluding Remarks

This paper has presented a sequential update algorithm for computing the stationary distribution vector in continuous-time upper block-Hessenberg Markov chains. The algorithm stops after finitely many iterations if Conditions 1 and 2 are satisfied. These conditions hold in any stable M/M/$s$ retrial queue and the BMAP/M/$\infty$ queues satisfying the mild condition (4.7). Furthermore, the algorithm would be applicable (under some mild conditions) to MAP/PH/$s$ retrial queues, BMAP/PH/$\infty$ queues, and their variants.

Appendix A Proof of Proposition 3.1

Let $\bm{T}_{n}^{*}$, $n\in\mathbb{Z}_{+}$, denote

[TABLE]

It then follows from (2.21), (3.9), and (A.1) that

[TABLE]

It also follows from (2.24), (A.1), and $\sum_{\ell=0}^{\infty}\bm{\pi}_{\ell}\bm{Q}_{\ell,n}=\bm{0}$ ($n\in\mathbb{Z}_{+}$) that

[TABLE]

Combining (A.2) and (A.3) yields (3.11). The proof has been completed.

Acknowledgments

The author thanks Mr. Masatoshi Kimura and Dr. Tetsuya Takine for providing the counterexample presented in Section 2.3. The author also thanks an anonymous referee for his/her valuable comments that helped to improve the paper.