The time-dependent expected reward and deviation matrix of a finite QBD process

TL;DR

This paper develops methods to compute the time-dependent expected reward and deviation matrix for finite QBD processes, aiding in performance analysis of queueing systems like MAP/PH/1/C.

Contribution

It introduces two approaches—matrix difference equations and perturbation theory—for calculating the deviation matrix of finite QBD processes, with practical numerical comparisons.

Findings

Both methods effectively compute the deviation matrix.

The recursive approach offers computational advantages.

Numerical examples demonstrate the methods' applicability.

Abstract

Deriving the time-dependent expected reward function associated with a continuous-time Markov chain involves the computation of its transient deviation matrix. In this paper we focus on the special case of a finite quasi-birth-and-death (QBD) process, motivated by the desire to compute the expected revenue lost in a MAP/PH/1/C queue. We use two different approaches in this context. The first is based on the solution of a finite system of matrix difference equations; it provides an expression for the blocks of the expected reward vector, the deviation matrix, and the mean first passage time matrix. The second approach, based on some results in the perturbation theory of Markov chains, leads to a recursive method to compute the full deviation matrix of a finite QBD process. We compare the two approaches using some numerical examples.

Equations

$$\mathbf{R}(0)$$

$${\mathbf{R}}^{\prime }(t)$$

$$\mathbf{R}(t)=(\mathbf{\pi }\mathbf{g})1t+D(t)\mathbf{g},$$

$$D(t)={\int }_0^t\left(e^{Qu}-1\mathbf{\pi }\right)du,$$

$$D={\int }_0^{\infty }\left(e^{Qu}-1\mathbf{\pi }\right)du,$$

$$Q{Q}^{\#}$$

$$W{Q}^{\#}$$

$$D=-Q^{\#}.$$

$$QD$$

$$\mathbf{\pi }D$$

$$D1=0.$$

$$V(t)={\int }_0^t{e}^{Qu}du=1\mathbf{\pi }t+D(t),$$

$$\mathbf{R}(t)=(\mathbf{\pi }\mathbf{g})1t+D(t)\mathbf{g},$$

$$\bar{\mathbf{R}}(t)=(\mathbf{\pi }\mathbf{g})1t+D\mathbf{g}.$$

$${\mathbf{R}}^{\prime }(t)=Q\mathbf{R}(t)+\mathbf{g}$$

$$\tilde{\mathbf{R}}(s)={\int }_0^{\infty }{e}^{-st}\mathbf{R}(t)dt,$$

$$\tilde{\mathbf{R}}(s)$$

$$=$$

$$\tilde{D}(s)=(1/s)(sI-Q{)}^{-1}-(1/s^2)1\mathbf{\pi }.$$

$$\tilde{D}(s)=(1/s)(sI-Q{)}^{-1}(I-1\mathbf{\pi }).$$

$$Q=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-5.0pt\left[\kern 0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{ \halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep& \kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&& \kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr 5.0pt\hfil\@arstrut$\scriptstyle\mbox{}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 1$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 2$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 3$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle C-1$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle C\\0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle B_{0}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle A_{1}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\\1$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle A_{-1}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle A_{0}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle A_{1}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\\2$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle A_{-1}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle A_{0}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle A_{1}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\\\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\\C-1$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle A_{-1}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle A_{0}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle A_{1}\\C$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle A_{-1}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle C_{0}\\$\hfil\kern 5.0pt\crcr}}}}\right]$}},$$

$$(B_0-sI){\tilde{\mathbf{R}}}_0(s)+A_1{\tilde{\mathbf{R}}}_1(s)$$

$$A_{-1}{\tilde{\mathbf{R}}}_{k-1}(s)+(A_0-sI){\tilde{\mathbf{R}}}_k(s)+A_1{\tilde{\mathbf{R}}}_{k+1}(s)$$

$$A_{-1}{\tilde{\mathbf{R}}}_{C-1}(s)+(C_0-sI){\tilde{\mathbf{R}}}_C(s)$$

$$A_{-1}+(A_0-sI)X+A_1{X}^2=0,$$

$$A_1+(A_0-sI)X+A_{-1}{X}^2=0.$$

$$\mathbb{E}\left[e^{s{\theta }_{k-1}}\text{mathds}1\{\phi ({\theta }_{k-1})=j\}|J(0)=k,\phi (0)=i\right]$$

$$\mathbb{E}\left[e^{s{\theta }_{k+1}}\text{mathds}1\{\phi ({\theta }_{k+1})=j\}|J(0)=k,\phi (0)=i\right].$$

$$H_0(s)=-(A_0-sI+A_{1}G(s)+A_{-1}\hat{G}(s){)}^{-1},$$

$$H_0=-(A_0+A_{1}G+A_{-1}\hat{G}{)}^{-1},$$

$$A_{-1}{\mathbf{u}}_{k-1}+A_0{\mathbf{u}}_k+A_1{\mathbf{u}}_{k+1}=-{\mathbf{g}}_k$$

$${\mathbf{u}}_k=G^k\mathbf{v}+{\hat{G}}^{C-k}\mathbf{w}+{\mathbf{\nu }}_k(C),$$

$${\mathbf{\nu }}_k(C)=\sum _{j=0}^{k-1}{G}^j{H}_0{\mathbf{g}}_{k-j}+\sum _{j=1}^{C-k}{\hat{G}}^j{H}_0{\mathbf{g}}_{k+j}.$$

Full text

We use two different approaches in this context. The first is based on the solution of a finite system of matrix difference equations; it provides an expression for the blocks of the expected reward vector, the deviation matrix, and the mean first passage time matrix. The second approach, based on some results in the perturbation theory of Markov chains, leads to a recursive method to compute the full deviation matrix of a finite QBD process. We compare the two approaches using some numerical examples.

Keywords: Finite quasi-birth-and-death process; expected reward; deviation matrix; matrix difference equations; perturbation theory.

1 Introduction

Our analysis of the expected reward in a finite QBD process is motivated by the problem of computing the expected amount of lost revenue in a MAP/PH/1/C queue over a finite time horizon $[0,t]$, given its initial occupancy. Since MAPs and PH distributions are the most general matrix extensions of Poisson processes and exponential distributions, respectively, we can think of this problem as a matrix generalisation of the similar analysis for the $M/M/C/C$ model considered in Chiera and Taylor [6] and the $M/M/1/C$ model in Braunsteins, Hautphenne and Taylor [4]. In these queueing models, customers are lost when they arrive to find $C$ customers already present. Assuming that each arriving customer brings a certain amount of revenue, we are interested in calculating the expected amount of revenue that the queue will lose over a finite time horizon $[0,t]$, as well as exploring the limit of the rate of losing revenue in the asymptotic regime.

Solving the expected lost revenue problem is important, for example, if we wish to find a way of managing a system where a number of, possibly different, queues share a number of servers. If it is feasible to reallocate servers from one queue to another every $t$ time units, then a rational method for performing the allocation is for each queue to observe its occupancy at time 0, calculate the expected revenue lost in time $[0,t]$ for a range of capacities, given its initial occupancy, and then to allocate the servers to minimise the total expected amount of lost revenue over $[0,t]$. At time $t$, the calculation can be performed again, based upon the occupancies at that time and a reallocation performed if it is optimal to do so.

The MAP/PH/1/C queue can be modelled as a finite QBD process with generator matrix $Q$, and levels $0,1,\ldots,C$ corresponding to the possible queue lengths. If $\boldsymbol{R}(t)$ is a vector containing the expected revenue lost in $[0,t]$ conditional on the initial state (level and phase) of the system, then computing $\boldsymbol{R}(t)$ reduces to solving a special case of the time-dependent version of Poisson’s equation of the form

[TABLE]

where $\boldsymbol{0}$ is a column vector of [math]’s, and $\boldsymbol{g}$ is a column vector containing the reward (loss) per unit of time in each state of the system. Since the MAP/PH/1/C system loses revenue only when it is at full capacity, the only non-zero entries of $\boldsymbol{g}$ in our motivating example are those corresponding to level $C$.

The solution of (1), given by

[TABLE]

where $\boldsymbol{\pi}$ is the stationary vector of the QBD process and $\boldsymbol{1}$ denotes the column vector of 1’s, involves the transient deviation matrix,

[TABLE]

see [4]. As $t\rightarrow\infty$, $D(t)$ converges to the deviation matrix $D$ discussed in Coolen-Schrijner and van Doorn [7], which corresponds to the group inverse of $-Q$, and the expected lost revenue function has a linear asymptote, $\boldsymbol{R}(t)\sim(\boldsymbol{\pi}\boldsymbol{g})\boldsymbol{1}\,t+D\boldsymbol{g}.$

After providing more detail on the reward function $\boldsymbol{R}(t)$ associated with a QBD process in Section 2, we tackle the computation of $\boldsymbol{R}(t)$ in transient and asymptotic regimes, and the corresponding matrices $D(t)$ and $D$, using two different approaches, each having some advantages in comparison to the other. In the first approach, developed in Section 3, we assume that the QBD process is non null-recurrent, and we place ourselves in the general context where the reward vector $\boldsymbol{g}$ is not restricted to any particular structure. We use systems of matrix difference equations to gain insight into the block-structure of the vector $\boldsymbol{R}(t)$ and of the matrices $D(t)$ and $D$. In addition, this method also provides us with the blocks of the matrix of the mean first passage times in the QBD process. We obtain simple expressions for the relevant quantities that highlight the role of the maximal capacity $C$. We also derive limiting results when the maximal capacity increases to infinity. This approach is effective if one wants to focus on particular blocks of $\boldsymbol{R}(t)$, $D(t)$ and $D$ rather than on the full matrices, as is the case in our motivating example. In practice, it also avoids dealing with large matrices: the size of the matrices involved in the expressions is at most twice the size of the phase space of the QBD process, regardless the value of $C$.

The second approach, described in Section 4, relies on some elegant results from the perturbation theory of Markov chains, and leads to a recursive formula expressing the full deviation matrices $D(t)$ and $D$ for a given value of $C$ in terms of the corresponding matrices for a system with capacity $C-1$.

The two approaches complement each other and their use may depend on the context: the first method is algebraic because we solve matrix difference equations, but the solutions have a probabilistic interpretation, while the second method is probabilistic but leads to a solution with an algebraic flavour.

In Section 5, we provide numerical illustrations, starting with our motivating MAP/PH/1/C example which loses revenue only at full capacity, and moving on to look at a case where the reward function $\boldsymbol{g}$ is non-zero for all levels of the system. We compare the computational complexity of the two approaches and show that there is always a threshold value of $C$ at which one method surpasses the other in terms of CPU time.

2 Background

Let $\{X(t):t\geq 0\}$ be an ergodic continuous-time Markov chain on a finite state-space, with generator $Q$ and stationary distribution $\boldsymbol{\pi}$. The deviation matrix of $\{X(t)\}$ is the matrix

[TABLE]

whose components may be written as ${D}_{ij}=\lim_{t\rightarrow\infty}[N_{ij}(t)-N_{\boldsymbol{\pi}j}(t)],$ where $N_{ij}(t)$ is the expected time spent in state $j$ during the interval of time $[0,t]$ given that the initial state is $i$, and $N_{\boldsymbol{\pi}j}(t)$ is the same quantity but conditional on the initial state having the distribution $\boldsymbol{\pi}$ (see Da Silva Soares and Latouche [8]).

The group inverse $A^{\#}$ of a matrix $A$, if it exists, is defined as the unique solution to $AA^{\#}A=A$, $A^{\#}AA^{\#}=A^{\#}$, and $A^{\#}A=AA^{\#}$. From Campbell and Meyer [5, Theorem 8.5.5], the group inverse $Q^{\#}$ of the infinitesimal generator $Q$ of any finite Markov chain is the unique solution of the two equations

[TABLE]

where $W=\lim_{t\rightarrow\infty}[\exp(Qt)]$. When it exists, the deviation matrix is related to the group inverse $Q^{\#}$ of $Q$ by the relation

[TABLE]

If $Q$ is irreducible, then $D$ is also the unique solution of the system

[TABLE]

see [7] for more detail. In addition, $D$ also satisfies

[TABLE]

We see that $D=\lim_{t\rightarrow\infty}D(t)$, where $D(t)$ is the transient deviation matrix defined by (3). Properties of the transient deviation matrix are discussed in [4].

Let $V_{ij}(t)$ be the expected cumulative time spent in state $j$ in the time interval $[0,t]$, given that the process starts in state $i$ at time [math], and define the matrix $V(t)=(V_{ij}(t))$. It follows that

[TABLE]

and $V(t)$ has the linear asymptote $\bar{V}(t)=\boldsymbol{1}\boldsymbol{\pi}\,t+D.$ If we associate a reward (or loss) $g_{j}$ per time unit when the Markov chain $\{X(t)\}$ occupies state $j$, and define the vector $\boldsymbol{g}:=(g_{j})$, then the expected cumulative reward up to time $t$, given that the chain starts in state $i$, is given by $R_{i}(t):=(V(t)\boldsymbol{g})_{i},$ so that the vector $\boldsymbol{R}(t):=(R_{i})$ satisfies

[TABLE]

and $\boldsymbol{R}(t)$ has the linear asymptote

[TABLE]

Observe that (10) is the solution of a finite horizon version of Poisson’s equation,

[TABLE]

with $\boldsymbol{R}(0)=\boldsymbol{0}$. In the Laplace transform domain where, for $\Re(s)>0$,

[TABLE]

(10) and (12) become respectively

[TABLE]

where $\tilde{D}(s)$ is the Laplace transform of the transient deviation matrix, given by

[TABLE]

Another equivalent expression, obtained after some algebraic manipulations, is given by

[TABLE]

In this paper, we are interested in computing the expected reward function $\boldsymbol{R}(t)$, and the associated deviation matrices $D(t)$ and $D$, for any $t\geq 0$, when the Markov chain corresponds to a finite (level-independent) QBD process. Such a process is a two-dimensional continuous-time Markov chain $\{\boldsymbol{X}(t)=(J(t),\varphi(t)),t\geq 0\}$, where the variable $J(t)$, taking values in $\{0,1,\ldots,C\}$, is called the level of the process at time $t$, and the variable $\varphi(t)$ taking values in $\{1,2,\ldots,n\}$ is called the phase of the process at time $t$. The generator of the QBD process has a block-tridiagonal form given by

[TABLE]

where the block matrices are of size $n\times n$.

The vectors $\boldsymbol{R}(t)$ and $\boldsymbol{g}$ corresponding to the QBD process are of length $(C+1)n$ and can be structured into block sub-vectors $\boldsymbol{R}_{k}(t)$ and $\boldsymbol{g}_{k}$ of length $n$, for $0\leq k\leq C$, corresponding to each level of the QBD process. The same block decomposition holds for the deviation matrices $D(t)$ and $D$, and we write $D_{k,\ell}(t)=\left(D_{(k,i)(\ell,j)}(t)\right)_{1\leq i,j\leq n}$ ($D_{k,\ell}=\left(D_{(k,i)(\ell,j)}\right)$ respectively) for the $(k,\ell)$th block of these matrices.

The next two sections address the following questions:

• (i)

How do we use the structure of the QBD process to compute the blocks of ${\boldsymbol{R}}(t)$, $D(t)$ and $D$?

• (ii)

How do the deviation matrices $D(t)$ and $D$ differ for two successive values of the capacity $C$?

3 A matrix difference equation approach

In this section, we answer Question (i). We deal with the transient regime first, followed by the asymptotic regime.

3.1 Transient regime

We analyse the transient behaviour of the expected reward function in the Laplace transform domain. An explicit expression for $\tilde{\boldsymbol{R}}(s)$ is given in (14) but this expression does not give any insight in the structure of the blocks $\tilde{\boldsymbol{R}}_{k}(s)$, which are conditional on the initial queue length, nor in the role of the maximum capacity $C$ in the solution. Instead, we rewrite (14) as a system of second-order matrix difference equations

[TABLE]

where $\boldsymbol{g}_{k}(s)=(1/s)\,\boldsymbol{g}_{k}$, and we focus on the computation of the blocks $\tilde{\boldsymbol{R}}_{k}(s)$, for $0\leq k\leq C$.

For any $s\geq 0$, let $G(s)$ be the minimal nonnegative solution to the matrix quadratic equation

[TABLE]

and let $\hat{G}(s)$ be the minimal nonnegative solution to

[TABLE]

The matrices $G(s)$ and $\hat{G}(s)$ can be computed numerically using any of the linear or quadratic algorithms discussed in Latouche and Ramaswami [12]. Assuming that the process is not restricted by an upper boundary, for $s>0$, the $(i,j)th$ entry of the matrix $G(s)$ contains the Laplace transform

[TABLE]

where $\theta_{k-1}=\inf_{t>0}\{J(t)=k-1\}$ and $\mathds{1}\{\cdot\}$ denotes the indicator function. Similarly, assuming that the process is not restricted by a lower boundary, the $(i,j)th$ entry of the matrix $\hat{G}(s)$ contains the Laplace transform

[TABLE]

Both $G(s)$ and $\hat{G}(s)$ are sub-stochastic, and therefore have a spectral radius strictly less than 1, as long as $s>0$. The first passage probability matrices $G$ and $\hat{G}$ correspond to $G(0)$ and $\hat{G}(0)$. When the unrestricted QBD process is transient, $G$ is sub-stochastic and $\hat{G}$ is stochastic, and conversely, when the unrestricted QBD process is positive recurrent, $G$ is stochastic and $\hat{G}$ is sub-stochastic. Finally, we let

[TABLE]

which is well defined for any $s>0$, and

[TABLE]

which is well defined in the non null-recurrent case.

In the sequel, we use the convention that an empty sum, such as $\sum_{j=1}^{0}$ or $\sum_{j=0}^{-1}$, is zero.

Lemma 3.1

For any $C\geq 1$ and $1\leq k\leq C$, the general solution of the second-order matrix difference equation

[TABLE]

is given by

[TABLE]

where $\boldsymbol{v}$ and $\boldsymbol{w}$ are arbitrary vectors, and

[TABLE]

Proof. First, using (20), (21) and (24) we can show that (26) is a solution to (25) for any arbitrary vectors $\boldsymbol{v},\boldsymbol{w}\in\mathbb{C}^{n}$. It remains to show that all solutions of (25) can be expressed in the form of (26).

Observe that there exists a nonsingular matrix $M$ such that

[TABLE]

with

[TABLE]

where $V$ is a non-singular square matrix of order $p$, and $W$ is a square matrix of order $q$ with $p+q=n$ and sp$(W)=0$. For instance, we may choose $M$ such that $J$ is the Jordan normal form of $\hat{G}$; in this case, $V$ contains all the blocks for the non-zero eigenvalues of $\hat{G}$ (which lie within the unit circle), and $W$ contain all the blocks for the zero eigenvalues. The matrices $M$ and $M^{-1}$ may be partitioned as

[TABLE]

where $L$ has dimension $n\times p$, $E$ has dimension $p\times n$, and $EL=I$, $FL=0$. Lemma 10 in Bini et al. [2] states that the general solution of (25) is given, for $k\geq 1$, by

[TABLE]

where $\boldsymbol{v}\in\mathbb{C}^{n}$ and $\boldsymbol{z}\in\mathbb{C}^{p}$ are arbitrary vectors, and

[TABLE]

with

[TABLE]

and $\nu$ the smallest integer such that $M^{\nu}=0$ ($\boldsymbol{\tau}_{k}=\boldsymbol{0}$ if $\hat{G}$ is invertible).

Since the eigenvalues of $V$ lie within the unit circle, the negative powers of $V$ appearing in (27) and (28) have entries that can take unbounded values, which can lead to numerical instabilities. Next, we rewrite (27) in a more convenient form in order to get rid of any negative power of $V$. First observe that (27) is equivalent to

[TABLE]

where the negative powers of $V$ appear in the term

[TABLE]

Letting $\boldsymbol{z}=\boldsymbol{y}+\sum_{j=1}^{C}V^{j}EH_{0}\,\boldsymbol{g}_{j}$, where $\boldsymbol{y}\in\mathbb{C}^{p}$ is an arbitrary vector, we obtain

[TABLE]

Let us fix a vector $\boldsymbol{y}\in\mathbb{C}^{p}$. Then, for any $k\leq C$, $LV^{-k}\,\boldsymbol{y}$ can be equivalently written as $\hat{G}^{C-k}\boldsymbol{w}$ with $\boldsymbol{w}=LV^{-C}\boldsymbol{y}$. Indeed,

[TABLE]

Therefore, we have

[TABLE]

and the general solution (27) takes the form

[TABLE]

where

[TABLE]

Finally, note that $\boldsymbol{g}_{\ell}$ is defined for $0\leq\ell\leq C$ only, so we can set $\boldsymbol{g}_{k+j}=\boldsymbol{0}$ for any $j>C-k$. With this, we have

[TABLE]

which shows that any solution to (25) can be written in the form (26). $\square$

The advantage of the solution (26) over the solution (27) from [2] is that it does not require any spectral decomposition of the matrix $\hat{G}$, nor any matrix inversion (as long as $C\geq k$). Since the spectral radii of $G$ and $\hat{G}$ are bounded by 1, all matrix powers involved in (26) are bounded, and the computation of the solution is therefore numerically stable.

Lemma 3.1 naturally extends to the Laplace transform domain, as stated in the next Corollary.

Corollary 3.2

For any $s>0$, $C\geq 1$, and $1\leq k\leq C$, the general solution of the second-order matrix difference equation

[TABLE]

is given by

[TABLE]

where $\boldsymbol{v}(s)$ and $\boldsymbol{w}(s)$ are arbitrary vectors, and

[TABLE]

$\square$

We use Corollary 3.2 to obtain a closed-form expression for $\tilde{\boldsymbol{R}}_{k}(s)$ in terms of the matrices $G(s)$ and $\hat{G}(s)$ of the QBD process.

Proposition 3.3

For any $s>0$, $C\geq 1$, and $0\leq k\leq C$, the Laplace transform of the expected reward function, conditional on the initial level $k$, is given by

[TABLE]

where $\boldsymbol{\nu}_{k}(s,C)$ is given by (31), and

[TABLE]

with

[TABLE]

In addition, $Z(s,C)$ can be written as

[TABLE]

where $\mathring{Q}(s,C)$ is the generator of the transient Markov chain obtained from the QBD process in which absorption can happen from any state at rate $s>0$, and restricted to levels 0 and $C$.

Proof. We apply Corollary 3.2 to (18) and obtain that the general solution of the second-order difference equation is given by (32) for $1\leq k\leq C-1$. We then specify the arbitrary vectors $\boldsymbol{v}(s,C)$ and $\boldsymbol{w}(s,C)$ using the boundary conditions (17) and (19). Injecting the general solution into (17) and (19) leads to the system of equations

[TABLE]

which can be rewritten as

[TABLE]

where $Z(s,C)$ is given by (33). To prove that the matrix $Z(s,C)$ is invertible, we now show that it can be written as the matrix product

[TABLE]

where $\mathring{Q}(s,C)$ is the generator of the transient Markov chain obtained from the QBD process in which absorption can happen from any state at rate $s>0$, restricted to levels 0 and $C$, and is therefore invertible, and the matrix inverse

[TABLE]

exists because sp$(G(s))<1$ and sp$(\hat{G}(s))<1$ for any $s>0$.

The generator $\mathring{Q}(s,C)$ of the transient Markov chain in which absorption can happen from any state at rate $s>0$, restricted to levels 0 and $C$, can be written as $\mathring{Q}(s,C)=\mathcal{A}(s)+\mathcal{B}(-\mathcal{C}(s))^{-1}\mathcal{D}$ with

[TABLE]

[TABLE]

We then have

[TABLE]

and to show that this is equal to $Z(s,C)$ given by (33), it suffices to show that

[TABLE]

Observe that

[TABLE]

where we have used (20) and (21). This directly leads to

[TABLE]

which, pre-multiplied by the matrix $\mathcal{B}$, provides (37). $\square$

Observe that the expression (32) for $\tilde{\boldsymbol{R}}_{k}(s)$ involves matrices of at most twice the size of the phase space (such as $Z(s,C)$). The function $\boldsymbol{R}_{k}(t)$ is obtained by taking the inverse Laplace transform of $\tilde{\boldsymbol{R}}_{k}(s)$, which can be done numerically using the method of Abate and Whitt [1] and the function $\mathsf{ilaplace.m}$ in Matlab.

The next corollary provides us with the limit as $C\rightarrow\infty$ of the result stated in Proposition 3.3. It is a direct consequence of the fact that the matrices $G(s)$ and $\hat{G}(s)$ are sub-stochastic for any $s>0$.

Corollary 3.4

Assume that the series $\sum_{j=1}^{\infty}\hat{G}(s)^{j}H_{0}(s)\,\boldsymbol{g}_{k+j}(s)$ converges for any $k\geq 0$. Then, for any $s>0$, and $k\geq 0$, the limit as $C\rightarrow\infty$ of the Laplace transform of the expected reward function, conditional on the initial level $k$, is given by

[TABLE]

where

[TABLE]

and

[TABLE]

$\square$

Note that we can write $\boldsymbol{\nu}_{0}(s,\infty)=\hat{G}(s)\boldsymbol{\nu}_{1}(s,\infty).$

In addition to the blocks $\tilde{\boldsymbol{R}}_{k}(s)$ of the expected reward function, Proposition 3.3 provides us with an explicit expression for the blocks $\tilde{D}_{k,\ell}(s)$ of the Laplace transform of the transient deviation matrix, as we show now.

Proposition 3.5

For $0\leq k\leq C$,

• •

if $\ell=0$,

[TABLE]

with

[TABLE]

• •

if $1\leq\ell\leq C-1$,

[TABLE]

with

[TABLE]

• •

if $\ell=C$,

[TABLE]

with

[TABLE]

where $Z(s,C)$ is given by (33).

Proof. We use (13) and Proposition 3.3 with, for $0\leq\ell\leq C$ and $1\leq i\leq n$, $\boldsymbol{g}_{\ell}=\boldsymbol{e}_{i}$ and $\boldsymbol{g}_{j}=\boldsymbol{0}$ for $j\neq\ell$. $\square$

The limit as $C\rightarrow\infty$ of Proposition 3.5 again follows easily, and provides us with an analytical expression for the blocks of the infinite matrix $\tilde{D}(s)$ corresponding to a QBD with no upper bound on the levels.

Corollary 3.6

For $k\geq 0$,

• •

if $\ell=0$,

[TABLE]

• •

if $\ell\geq 1$,

[TABLE]

with

[TABLE]

$\square$

3.2 Asymptotic regime

Next, we concentrate on the asymptotic properties of the expected reward function ${\boldsymbol{R}}_{k}(t)$ for large values of $t$. By decomposing (11) into blocks, we have

[TABLE]

We determine now an explicit expression for the blocks $D_{k,\ell}$ of the deviation matrix for $0\leq k,\ell\leq C$, together with the mean first passage times to any level $\ell$ in the QBD process. Our method is based on the relationship between the entries of the deviation matrix and mean first passage times, and involves the solution of finite systems of matrix difference equations similar to the ones we solved in Section 3.1. Indeed, we can write

[TABLE]

where $M_{(k,i)(\ell,j)}$ is the mean first entrance time to $(\ell,j)$ from $(k,i)$, and $M_{\boldsymbol{\pi}(\ell,j)}$ is the mean first entrance time to $(\ell,j)$ if the state at time 0 has the stationary distribution $\boldsymbol{\pi}$. Define the block-matrices $M_{k,\ell}=\left(M_{(k,i)(\ell,j)}\right)_{1\leq i,j\leq n}$, for $0\leq k,\ell\leq C$. In matrix form, we have

[TABLE]

Let us fix level $\ell$ and define $\boldsymbol{m}^{(j)}_{k,\ell}=M_{k,\ell}\,\boldsymbol{e}_{j}$, the $j$th column of $M_{k,\ell}$, for $1\leq j\leq n$ and $0\leq k\leq C$. So $(\boldsymbol{m}^{(j)}_{k,\ell})_{i}=M_{(k,i)(\ell,j)}$. Since the diagonal of $M_{\ell,\ell}$ is null we must have $(\boldsymbol{m}^{(j)}_{\ell,\ell})_{j}=0$. We introduce the notation $\bar{A}_{-1}^{(j)}$ and $\bar{A}_{1}^{(j)}$ for the matrices obtained by replacing the $j$th row in $A_{-1}$ and $A_{1}$, respectively, by $\boldsymbol{0}^{\top}$, and $\bar{B}_{0}^{(j)}$, $\bar{A}_{0}^{(j)}$, $\bar{C}_{0}^{(j)}$ for the matrices obtained by replacing the $j$th row in $B_{0}$, $A_{0}$ and $C_{0}$, respectively, by $-\boldsymbol{e}_{j}^{\top}$.

Proposition 3.7

For any fixed level $0\leq\ell\leq C$ and any phase $1\leq j\leq n$, the vectors $\boldsymbol{m}^{(j)}_{k,\ell}$ satisfy the system of matrix second-order difference equations

[TABLE]

for $1\leq k\leq\ell-1$ and $\ell+1\leq k\leq C-1$, with the following boundary conditions, depending on the value of $\ell$:

• •

for $\ell=0$,

[TABLE]

• •

for $1\leq\ell\leq C-1$,

[TABLE]

• •

for $\ell=C$,

[TABLE]

Proof. For a fixed value of $\ell$, let $m^{(j)}_{(k,i)}=(\boldsymbol{m}^{(j)}_{k,\ell})_{i}$, for $1\leq i\leq n$. By conditioning on the epoch where the process first leaves the state $(k,i)$, we obtain for $1\leq k\leq C-1$, $\ell+1\leq k\leq C-1$, and $1\leq j\leq n$,

[TABLE]

which in matrix form gives (41).

A similar argument leads to the boundary equations (43), (44), (46), and (47). Finally, the boundary equations (42), (45), and (48) are obtained by adding the constraint that when $k=\ell$, $m^{(j)}_{(\ell,j)}=0$. $\square$

For $0\leq k\leq C$, define the vectors

[TABLE]

The next theorem provides an explicit expression for the columns of the mean first passage time matrices $M_{k,\ell}$.

Proposition 3.8

In the non null-recurrent case, for any $0\leq k\leq C$, the vector $\boldsymbol{m}^{(j)}_{k,\ell}$ has the following explicit expression, depending on the value of $\ell$:

• •

if $\ell=0$,

[TABLE]

where $\boldsymbol{\mu}_{k}(C)$ is given by (50) and

[TABLE]

with

[TABLE]

where $\mathring{Q}^{(j)}(C)$ is the generator of the transient Markov chain obtained from the QBD process in which absorption can happen from state $(0,j)$ at rate $1$, and restricted to levels 0 and $C$.

• •

if $1\leq\ell\leq C-2$,

[TABLE]

where $\boldsymbol{\mu}_{k}(C)$ is given by (50) and

[TABLE]

with

[TABLE]

where $\mathring{Q}^{(j)}(\ell,C)$ is the generator of the transient Markov chain obtained from the QBD process in which absorption can happen from state $(\ell,j)$ at rate $1$, and restricted to levels 0, $\ell$, $\ell+1$, and $C$.

• •

if $\ell=C-1$,

[TABLE]

where $\boldsymbol{\mu}_{k}(C)$ is given by (50) and

[TABLE]

with

[TABLE]

where $\mathring{Q}^{(j)}(C)$ is the generator of the transient Markov chain obtained from the QBD process in which absorption can happen from state $(C-1,j)$ at rate $1$, and restricted to levels 0, $C-1$ and $C$.

• •

if $\ell=C$,

[TABLE]

where $\boldsymbol{\mu}_{k}(C)$ is given by (50) and

[TABLE]

with

[TABLE]

where $\mathring{Q}^{(j)}(C)$ is the generator of the transient Markov chain obtained from the QBD process in which absorption can happen from state $(C,j)$ at rate $1$, and restricted to levels 0 and $C$.

Proof. The result follows from Lemma 3.1 applied to the system of difference equations (41). The arbitrary vectors are then determined using the boundary conditions given in Proposition 3.7.

When $\ell=0$ and $\ell=C$, there are two boundary equations, determining two arbitrary vectors.

When $1\leq\ell\leq C-2$, the solution depends on whether $0\leq k\leq\ell$ or $\ell+1\leq k\leq C$. There are four boundary conditions, namely the three equations described in Proposition 3.7, in addition to one boundary equation obtained by taking $k=\ell+1$ in (41). This determines the four arbitrary vectors.

When $\ell=C-1$ the solution depends on whether $0\leq k\leq C-1$ or $k=C$. There are three boundary equations, as described in Proposition 3.7, which determine the three arbitrary vectors.

The decomposition of the matrices $Z^{(j)}(C)$ and $Z^{(j)}(\ell,C)$ into the product of a non-conservative generator and the matrices involving $G$ and $\hat{G}$ follows from the same arguments as those used in the proof of Proposition 3.3. Absorption from state $(\ell,j)$ at rate 1 comes from the definition of $\bar{B}_{0}^{(j)}$, $\bar{A}_{0}^{(j)}$, and $\bar{C}_{0}^{(j)}$. The matrices involving $G$ and $\hat{G}$ are invertible because sp$(G)<1$ or sp$(\hat{G})<1$ in the non null-recurrent case. $\square$

Observe that in the transient case the series $\sum_{j=1}^{\infty}\hat{G}^{j}H_{0}\,\boldsymbol{1}$ diverges, and the limit as $C\rightarrow\infty$ of $\boldsymbol{\mu}_{k}(C)$ given in (50) is infinite, while in the positive recurrent case,

[TABLE]

This leads to the following corollary to Proposition 3.8 for the mean first passage times in an unrestricted positive recurrent QBD process:

Corollary 3.9

In the positive recurrent case, for any $k\geq 0$, depending on the value of $\ell$, the vector $\boldsymbol{m}^{(j)}_{k,\ell}$ has the explicit expression:

• •

if $\ell=0$,

[TABLE]

where $\boldsymbol{\mu}_{k}(\infty)$ is given by (59), and

[TABLE]

• •

if $\ell\geq 1$,

[TABLE]

where $\boldsymbol{\mu}_{k}(\infty)$ is given by (59) and

[TABLE]

with

[TABLE]

$\square$

To complete the characterisation of the block matrices $D_{k,\ell}$ using (40), it remains for us to compute the blocks $\boldsymbol{\pi}_{x}$ of the stationary distribution of the QBD process for $0\leq x\leq C$. This can be done following Theorem 10.3.2 in [12] or Hajek [9], adapted to the continuous-time setting. This involves the matrices $R$ and $\hat{R}$ of rates of sojourn in level $k+1$, respectively $k-1$, per unit of the local time in level $k$ in the corresponding unrestricted QBD process. These matrices can be expressed in terms of $G$ and $\hat{G}$ as $R=A_{1}(-(A_{0}+A_{1}G))^{-1}$ and $\hat{R}=A_{-1}(-(A_{0}+A_{-1}\hat{G}))^{-1}$, see [12].

Corollary 3.10

In the non null-recurrent case, the stationary distribution of the QBD process with transition matrix (16) is given by

[TABLE]

where $(\boldsymbol{v}_{0},\boldsymbol{v}_{C})$ is the solution of the system

[TABLE]

and

[TABLE]

4 A perturbation theory approach

In this section, we answer Question (ii) on Page (ii) and we provide an expression for the deviation matrix of a QBD with maximal level $C$ in terms of the deviation matrix of the same QBD with capacity $C-1$, in transient and asymptotic regimes. To highlight the number of levels in the QBDs, we shall denote by $Q^{(C)}$ the generator (16) of a QBD with maximum level $C$, and by $D^{(C)}(t)$, $\tilde{D}^{(C)}(s)$, and $D^{(C)}$ the corresponding deviation matrices, for $C\geq 1$.

The key idea behind the derivation of the recursions is the decomposition of the generator $Q^{(C)}$ into a sum of another generator and a perturbation matrix, so that $Q^{(C)}$ can be seen as a block-element updating problem, the block version of the element-updating problem described in Langville and Meyer [11]. Precisely, we write

[TABLE]

where

[TABLE]

where $E_{C-1}^{(C)}$ has dimensions $n(C+1)\times n$ and $\Delta^{(C)}$ has dimensions $n\times n(C+1)$. Observe that ${T}^{(C)}$ is the generator of a reducible transient structured Markov chain, and its stationary distribution is given by

[TABLE]

where $\boldsymbol{\pi}^{(C-1)}$ is the stationary distribution of $Q^{(C-1)}$. The group inverse of $T^{(C)}$ exists by [5, Theorem 8.5.5], and may be expressed as a function of $\boldsymbol{\pi}^{(C-1)}$ and the deviation matrix $D^{(C-1)}=-(Q^{(C-1)})^{\#}$, as shown in the next lemma.

Lemma 4.1

The group inverse of $T^{(C)}$ is given by

[TABLE]

where $M^{(C)}$ is an $n\times nC$ matrix defined as

[TABLE]

Proof. We check that $(T^{(C)})^{\#}$ is the group inverse of $T^{(C)}$ by direct verification of (5, 6). $\square$

First, we use the decomposition (61) to derive an expression for $\boldsymbol{\pi}^{(C)}$ in terms of $\boldsymbol{\pi}^{(C-1)}$.

Lemma 4.2

The stationary distribution of $Q^{(C)}$ may be expressed in terms of the stationary distribution of $Q^{(C-1)}$ as

[TABLE]

Proof. Our argument is a slight modification of the proof of Rising [14, Lemma 5.1], which deals with irreducible Markov matrices. Let $A$ and $\widetilde{A}$ be two generators on the same state space, assume that each has a single irreducible class of states, not necessarily the same, so that each has a unique stationary probability vector. Denote the stationary distribution vectors as $\boldsymbol{\alpha}$ and $\widetilde{\boldsymbol{\alpha}}$, respectively.

By (7), $\widetilde{\boldsymbol{\alpha}}AA^{\#}=\widetilde{\boldsymbol{\alpha}}-\boldsymbol{\alpha}$ and so

[TABLE]

since $\widetilde{\boldsymbol{\alpha}}\widetilde{A}=\boldsymbol{0}$. If $A$ and $\widetilde{A}$ are irreducible, it results from [14, Lemma 5.1] that $I+(\widetilde{A}-A)A^{\#}$ is non-singular. A key argument in [14] is that the kernel of both $A$ and $\widetilde{A}$ is limited to Span($\{\boldsymbol{1}\}$), and so we repeat it verbatim to conclude that $I+(\widetilde{A}-A)A^{\#}$ is non-singular in our case as well. Thus,

[TABLE]

and the lemma is proved, once we set $A=T^{(C)}$ and $\widetilde{A}=Q^{(C)}$. $\square$

We thus see that $\boldsymbol{\pi}^{(C)}$ can be expressed in terms of $\boldsymbol{\pi}^{(C-1)}$ and $D^{(C-1)}$ through $(T^{(C)})^{\#}$. In the next two sections, we use the decomposition (61) to derive recursive formulae for $\tilde{D}^{(C)}(s)$ and $D^{(C)}$.

4.1 Time-dependent deviation matrix

Recall that the Laplace transform of the transient deviation matrix of $Q^{(C)}$ is given by

[TABLE]

We can express $\tilde{D}^{(C)}(s)$ in terms of $\tilde{D}^{(C-1)}(s)$ thanks to the fact that $(sI-Q^{(C)})^{-1}$ can be expressed in terms of $(sI-Q^{(C-1)})^{-1}$ via the Sherman-Morrison-Woodbury (SMW) formula, which we restate here.

Let $M$ be a finite non-singular matrix of order $m$ and let $U$ and $V$ be matrices of dimension $m\times k$ and $k\times m$, respectively. If $M+UV^{\top}$ is non-singular, then $I+V^{\top}M^{-1}U$ is non-singular as well and

[TABLE]

As a direct consequence of this formula, we have the following

Corollary 4.3

The inverse of $(sI-Q^{(C)})$ satisfies the recursive equation

[TABLE]

where $\Delta^{(C)}$ is given in (70), and

[TABLE]

where $M^{(C)}$ is given in (72).

Proof. We use the decomposition (61) and apply (75) with $M=sI-T^{(C)}$, $U=E_{C-1}^{(C)}$, and $V=-\Delta^{(C)\top}$. The expression for $(sI-T^{(C)})^{-1}$ can be checked by verifying that $(sI-T^{(C)})^{-1}(sI-T^{(C)})=I$ using (66). $\square$

The inverse $(sI-Q^{(C)})^{-1}$ is expressed in terms of $(sI-Q^{(C-1)})^{-1}$ and, by (74),

[TABLE]

Furthermore, since $\boldsymbol{\pi}^{(C)}$ can be expressed in terms of $\boldsymbol{\pi}^{(C-1)}$ and $D^{(C-1)}$ by Lemmas 4.1 and 4.2, we conclude that $\tilde{D}^{(C)}(s)$ may be expressed in terms of $\tilde{D}^{(C-1)}(s)$, $\boldsymbol{\pi}^{(C-1)}$, and $D^{(C-1)}$. The recursive computation of the transient deviation matrices $\tilde{D}^{(C)}(s)$ may therefore be done together with the recursive computation of the stationary distribution vectors $\boldsymbol{\pi}^{(C)}$ and the deviation matrices ${D}^{(C)}$, for $C\geq 1$.

4.2 Deviation matrix

Now, we focus on obtaining an expression for $D^{(C)}$ in terms of $D^{(C-1)}$ using the decomposition (61). We start from Theorem 3.1 in [11], which gives an expression for the group inverse of a one-element updated Markov chain generator, and we extend it to a one-block update. Starting from a generator $Q$, we consider a new generator $\tilde{Q}=Q+A$ where $A$ is a matrix with one non-zero block-row only. We denote by E_{K}^{(C)}=[\begin{array}[]{cccccc}0&\cdots&0&I&\cdots&0\end{array}]^{\top} the block-column with identity as $K$-th block and $C$ zero blocks elsewhere, and by $P$ a block-row vector containing $C+1$ blocks, then $A=E_{K}^{(C)}P$ for some $K$.

Proposition 4.4

Let $Q$ be the generator of a Markov process with deviation matrix $D$ and with stationary distribution $\boldsymbol{q}$. Suppose that $\tilde{Q}=Q+E_{K}^{(C)}P$ is the generator of an irreducible Markov process. The deviation matrix $\tilde{D}$ of $\tilde{Q}$ is given by

[TABLE]

where $\tilde{\boldsymbol{\pi}}$ is the stationary distribution of $\tilde{Q}$.

Proof. Observe first that the inverse in (77) is well defined by [14, Lemma 5.1]. We first show (77) by direct verification of (7, 8). Replacing $\tilde{D}$ with the right-hand side of (77), we obtain

[TABLE]

since $\boldsymbol{q}(I-E_{K}^{(C)}PD)^{-1}=\tilde{\boldsymbol{\pi}},$ by (73), where $A$ is replaced by $Q$ and $\widetilde{A}$ by $Q+E_{K}^{(C)}P$. Next,

[TABLE]

by the SMW formula (75), and injecting (79) into (77) provides (78). $\square$

Remark 4.5

Equation (78) can be seen as an extension of the SMW formula (75) to singular matrices.

Remark 4.6

As observed in [11], the cost of straightforward computation of (77) may be high if one updates more than one row, as it requires the inversion of the matrix $I-E_{K}^{(C)}PD$ whose size is the same as $\tilde{Q}$. Equation (78) reduces the computational cost because the size of $I-PDE_{K}^{(C)}$ is just one block.

Since the decomposition of $Q^{(C)}$ given in (61) is a block-element-updating of the matrix $T^{(C)}$ by $E_{C-1}^{(C)}\Delta^{(C)}$, the recursive formula for the deviation matrix $D^{(C)}$ of $Q^{(C)}$ follows from Proposition 4.4.

Proposition 4.7

The deviation matrix $D^{(1)}$ of the QBD process with generator

[TABLE]

is given by

[TABLE]

For $C\geq 2$, the deviation matrix $D^{(C)}$ of the QBD process with generator $Q^{(C)}$ is recursively given by

[TABLE]

where $\Delta^{(C)}$ is defined in (70), and where $(T^{(C)})^{\#}$ is given in terms of $D^{(C-1)}$ in (71).

Proof. The group inverse of $Q^{(1)}$ comes directly from the relationship between the group inverse of a generator and Kemeny and Snell’s fundamental matrix [10], given in Theorem 3.1 in [13].

For any $C\geq 2$, Equation (82) corresponds to Equation (78) in Proposition 4.4 with $\tilde{D}=D^{(C)}$, $\tilde{\boldsymbol{\pi}}=\boldsymbol{\pi}^{(C)}$, $D=-(T^{(C)})^{\#}$, $K=C-1$, and $P=\Delta^{(C)}$. $\square$

5 Numerical illustrations

In this section, we compute the expected loss revenue function for various MAP/PH/1/C queues using the results of Section 3. We then compare the performance of the two different approaches developed in Sections 3 and 4 by computing the deviation matrix of arbitrary finite QBD processes with different numbers of phases and levels.

5.1 The MAP/PH/1/C queue.

Recall that our motivation behind the computation of the expected reward function and the deviation matrices of a finite QBD process stems from our desire to compute of the expected amount of revenue lost in a MAP/PH/1/C queue. In such a process,

• (i)

the arrival MAP is characterised by the matrices $D_{0}$ and $D_{1}$ with $n_{1}$ phases;

• (ii)

the service time is distributed according to a PH$(\boldsymbol{\tau},T)$ distribution of order $n_{2}$, with $\boldsymbol{t}=-T\boldsymbol{1}$.

The MAP/PH/1/C system is a finite QBD process $\{\boldsymbol{X}(t)=(J(t),\boldsymbol{\varphi}(t)),t\geq 0\}$ where

• •

$0\leq J(t)\leq C$ represents the number of customers in the system at time $t$,

• •

$\boldsymbol{\varphi}(t)=(\varphi_{1}(t),\varphi_{2}(t))$ where $0\leq\varphi_{1}(t)\leq n_{1}$ is the phase of the MAP at time $t$, and $0\leq\varphi_{2}(t)\leq n_{2}$ is the phase of the PH distribution at time $t$.

The generator of that QBD process has a block-tridiagonal form given by (16) where the block matrices are of size $n=n_{1}\,n_{2}$ and are given by

[TABLE]

and

[TABLE]

where $\otimes$ and $\oplus$ denote the Kronecker product and sum, respectively.

We consider a simple example where the MAP of arrivals corresponds to a PH renewal process with

[TABLE]

and the service time is PH$(\boldsymbol{\tau},T)$ with

[TABLE]

Let $A=A_{-1}+A_{0}+A_{1}$ be the phase transition matrix associated with the QBD process, and let $\boldsymbol{\alpha}$ be the stationary vector of $A$. Since in this example $\boldsymbol{\alpha}A_{-1}\boldsymbol{1}<\boldsymbol{\alpha}A_{1}\boldsymbol{1}$, we are in a high-blocking system. This would correspond to a transient QBD process if there were no bound on the number of levels [12].

We first assume that each customer accepted in the system generates $\theta$ units of revenue. The system loses revenue when it is in level $C$ and a customer arrives, who is rejected, and does not lose any revenue otherwise. The expected reward function $\boldsymbol{R}(t)$ records the total expected amount of lost revenue over the finite time horizon $[0,t]$. In this particular case, the reward vector $\boldsymbol{g}$ records the loss per unit time, per state, and we have $\boldsymbol{g}_{k}=\boldsymbol{0}$ for $0\leq k<C$ and $\boldsymbol{g}_{C}=\theta A_{1}\boldsymbol{1}$. Results when $C=5$ and $\theta=1$ are shown in the upper panel of Figure 1 where we have assumed that the initial phase follows the distribution $\boldsymbol{\alpha}$.

Swapping the PH distribution of the arrival process and service time leads to a low-blocking system, whose results are shown in the lower panel of Figure 1. We clearly see that the expected lost revenue is much lower in the low-blocking system than in the high-blocking system, and in both cases it increases with the initial queue size, as expected.

We now consider a case where each customer brings a fixed amount of revenue $\gamma$ per time unit when in the system, in addition to the fixed $\theta$ units of revenue when entering the system. The reward vector then becomes $\boldsymbol{g}_{k}=\theta A_{1}\boldsymbol{1}+\gamma k\boldsymbol{1}$ for $0\leq k\leq C-1$ and $\boldsymbol{g}_{C}=\gamma C\boldsymbol{1}$, and we are now interested in the expected amount of revenue gained in $[0,t]$. The results are shown in Figure 2 for a high-blocking system (upper panel) and a low-blocking system (lower panel).

5.2 Comparison of the two approaches.

The matrix difference approach developed in Section 3 allows us to obtain the blocks of the expected reward vector and of the deviation matrices, with the advantage that the matrices involved in the expressions are at most twice the size of the phase space. However, that method requires the numerical computation of the matrices $G(s)$ and $\hat{G}(s)$ for $s\geq 0$; in addition, if we wish to compute the blocks of the deviation matrix, we first need to compute the blocks of the mean first passage time matrix, column by column. In contrast, the perturbation theory approach developed in Section 4 allows us to obtain recursively the whole deviation matrix, but at the cost of dealing with operations on larger matrices.

We measured the numerical efficiency of the two approaches by comparing the CPU time corresponding to the computation of the last block-column of the deviation matrix for arbitrary QBD processes with a different number of phases and a different maximal level $C$. To that aim, we generated QBD processes with random entries, with $n=2,\ldots,5$ and $C=1,\ldots,100$, and for each process, we ran the two methods 100 times and took the average CPU time. We repeated the same simulation several times and observed similar behaviour. The results of one experiment are shown in Figure 3. The CPU time increased linearly with $C$ for the matrix difference approach, while the growth is exponential for the perturbation theory approach. We see that there is always a threshold value $C^{*}$ such that the perturbation theoretical approach outperforms the matrix difference approach when $C<C^{*}$, and conversely when $C>C^{*}$. Note that the threshold value does not seem to show a monotonic behaviour when the number of phases $n$ increases.

In conclusion, if the objective is to compute parts of the deviation matrix, the perturbation theory approach is preferable for relatively small values of the maximal capacity, while the matrix difference approach is clearly more efficient for large values of $C$.

Acknowledgements

Sophie Hautphenne and Peter Taylor would like to thank the Australian Research Council (ARC) for supporting this research through Discovery Early Career Researcher Award DE150101044 and Laureate Fellowship FL130100039, respectively. In addition, they both acknowledge the support of the ARC Center of Excellence for Mathematical and Statistical Frontiers (ACEMS). \ Sarah Dendievel and Guy Latouche thank the Ministère de la Communauté française de Belgique for funding this research through the ARC grant AUWB-08/13-ULB 5. Sarah Dendievel would also like to thank the Methusalem program of the Flemish Community of Belgium for supporting this work in part.