Analyzing a single hyper-exponential working vacation queue from its governing difference equation

Abstract

As the queue becomes exhausted, different maintenance tasks can be performed according to the fatigue load and wear degree of the service equipment. At the same time, considering the customer's sensitivity to time delay, the service facility will not completely remain inactive during the maintenance period. To describe this objectively existing phenomenon arising in the waiting line system, we consider a hyper-exponential working vacation queue with a batch renewal arrival process. Through the calculation of the well-structured roots of the associated characteristic equation, the shift operator method in the theory of difference equations and the supplementary variable technique for stochastic modeling plays a central role in the queue-length distribution analysis. Comparison with other ways to analyze queueing models, the advantage of our approach is that we can avoid deriving the complex transition probability matrix of the queue-length process embedded at input points. The feasibility of this approach is verified by extensive numerical examples.

1. Introduction

The queueing system with server vacations is useful to model a system in which the server has an additional task during its idle period. The additional task may represent the server's working on some supplementary jobs, performing service equipment maintenance inspections and repairs, or server's rest after queue exhaustion. Since such a queue has broad applicability in analyzing the performance of industrial production systems and data communication networks, it has attracted many researchers’ attention over the past several decades. Various vacation policies provide more flexibility for optimal design and operating control of the systems. The working vacation policy introduced by Servi and Finn [23] is a kind of semi-vacation policy. It is characterized by the feature that the service facility works at a lower service rate rather than completely stopping service during the vacation period. At first glance, the concept of working vacation is a bit of an oxymoron because work and vacation are essentially two different things. But when we combine them into a waiting line system arising from the manufacturing environment, it might be a perfect way to reduce customer service response time because we can maintain the production equipment during the vacation period without suspending production completely. In the last twenty years, numerous researchers, including Wu and Takagi [26], Liu et al. [19], Li et al. [18], Zhang and Hou [28,29] , Selvaraju and Goswami [22], Gao and Yao [11], Lee and Kim [16], Ma et al. [20], used different methods to study this kind of queue under the assumption that customers arrive at a service facility according to a Poisson stream. Since our research is mainly concerned with the general renewal input working vacation queue, we only give a brief literature review in this area and point out the limitations of the current study to clarify our work's motivation.

So far, the matrix-geometric approach and the embedded Markov Chain technique are the mainstream methods to study the GI/M/1 type queue, and there is no exception for the working vacation queue. Baba [1] applied the matrix-geometric method to investigate a GI/M/1 queue with multiple exponential working vacations. Li et al. [17] used the same way to study the discrete-time version of the above model. Furthermore, both in continuous and discrete-time cases, Chae et al. [7] discussed general input queues with a single working vacation using the embedded Markov chain. Considering the fact that customers arrive in batches instead of individually, Guha and Banik [14] and Guha et al. [15] also analyzed the renewal input batch arrival queue under single and multiple exponential working vacation policies. Recently, Yu [27] extended the above work so that the bulk service rule is included. After carefully reading these papers mentioned above, we may note that each one contains a more complex transition probability analysis for the queue-length process embedded at the pre-arrival epoch. The corresponding probabilistic arguments involved in these studies are cumbersome to implement and error-prone. Additionally, for the general input working vacation queue, except for the work done by Chen et al. [9], nearly all other authors assume that the server's vacation time is exponentially distributed. We think that this is not a very reasonable assumption for many cases in the industrial environment. Under such a situation, when the waiting line system becomes exhausted, we usually decide what maintenance needs to be done according to the fatigue state of service equipment. Maintenance may represent the actions of replacement, repair, adjustment, overhaul, improvement, and checking of the service equipment. Since each action corresponds to a degree of complexity, the maintenance time cannot be assumed to follow a common exponential distribution. Thus, to give a more realistic model assumption, we intend to extend the exponential vacation time to hyper-exponential vacation time in this research. Moreover, to provide an alternative yet simple problem-solving methodology, we try to use an approach based on the theory of difference equations to carry out a comprehensive analysis of GI$^X$/M/1 single working vacation queue. Here, we must admit that the recent work done by Barbhuiya and Gupta [2–5] gives us some basic ideas to complete the model analysis. We will see that the method adopted in this paper can make us get rid of the chain matrix analysis. So for most people, this is an easy to accept method.

The rest of this paper is organized as follows. In Section 2, we describe the mathematical queueing model in detail. A set of differential-difference equations that represent the dynamics of the queue-length process is developed in Section 3. In Section 4, queue-length probabilities at different epochs are derived explicitly by solving simultaneous nonhomogeneous difference equations. Section 5 is devoted to the sojourn time of a randomly selected customer in an arriving batch. A three-step algorithm for computing the queue-length distribution is summarized in Section 6. To validate our computational algorithm, we further provide several typical numerical examples in this section. Finally, conclusions and future scopes are presented in Section 7.

2. Model formulation

The model is defined by making the following assumptions.

1. (1) Consider a single-station queueing system where customers arrive in batches according to a renewal process with independent identically distributed (i.i.d.) inter-batch arrival times having a common cumulative distribution function $A(t)$, and probability density function $a(t)$. Let $a^*(s)=\int _0^\infty e^{-st}\,{\rm d}A(t)$ be the Laplace–Stieltjes transform (L.S.T.) of $A(t)$ and let the mean inter-batch arrival time be denoted by $1/\lambda$. Differentiation of the L.S.T. with respect to $s$ is justified and yields $1/\lambda =-({{\rm d}}/{{\rm d}s})a^*(s)|_{s=0}<\infty$.

2. (2) At every arrival epoch, a batch of $k$ customers arrives with probability $g_k$. For mathematical convenience and from a more realistic point of view, we assume that the maximum batch size is equal to $b$. Consequently, the probability generating function of the sequence $\{g_k,k=1,2,\ldots,b\}$ is $G(z)=\sum _{k=1}^b g_kz^k$ with the first moment $\overline {g}=\sum _{k=1}^b kg_k$.

3. (3) In a normal busy period, the service times are i.i.d. exponential random variables with mean $1/\mu _0$. When the system becomes empty, the server goes on vacation for a random duration, where the vacation time has an $h$-phase hyper-exponential distribution, which can be properly represented by a probabilistic mixture of exponential distributions. The distribution function of the vacation time $V$ is defined as $\Pr \{V\le t\}=V(t)=\sum _{j=1}^h\alpha _j(1-e^{-\theta _jt})$. In other words, at the end of a service, if no customer is left in the system, the server takes type $j$ vacation with probability $\alpha _j$, then the vacation time distribution is ${\rm Exp}(\theta _j)$, $j=1,2,\ldots,h$.

4. (4) Customers arriving during a type $j$ vacation will be served at a lower service rate, where the service time obeys the exponential distribution with rate $\mu _j$ ($\mu _j<\mu _0$, $j=1,2,\ldots,h$). Upon completion of service at a lower rate, if the vacation is not over, and no customer is waiting for service, the server will continue the current vacation. On the contrary, if there is at least one customer in the system, the server will keep this service mode.

5. (5) As this type of $j$ vacation gets over, the server turns to the normal working mode immediately. The service of a customer being served will be interrupted and restarted from the beginning in the normal busy period. Alternatively, if no customers are found in the queue at the end of a vacation, the server remains idle and is ready to serve new arrivals at a normal service rate $\mu _0$.

6. (6) We further assume that the inter-batch arrival times, service times, and vacation times are mutually independent. A necessary and sufficient condition for model stability is $\rho =\lambda \overline {g}/\mu _0<1$ (see the proofs of Lemma 1 and the analysis of Eq. (35)).

3. Governing difference equation

Queueing systems with general inter-batch arrival time distribution and hyper-exponential vacation time are difficult to analyze mathematically due to the queueing process being non-Markovian. To enable the system to be characterized as a Markov system, the following random variables will be used for the development of our model. Let $N(t)$ and $\xi (t)$ denote the number of customers in the system (including the one in service) and the state of the server at time $t$, respectively. Here

$$\xi(t)=\left\{\begin{array}{ll} 0, & \text{if the server is in normal busy period at time}\ t,\\ j, & \text{if the server is on type} \ j\ \text{vacation at time}\ t,\ j=1,2,\ldots,h. \end{array}\right.$$

With the inclusion of the supplementary variable technique, and together with the remaining inter-batch arrival time $\tilde {A}(t)$ at time $t$, we may obtain a tri-variate Markov process $\{N(t),\xi (t),\tilde {A}(t)\}$. Furthermore, to establish the dynamic model of the above Markov process, let us define some probabilities as below

\begin{align*} & P_i(x,t)\,{\rm d}x=\Pr\{N(t)=i,\xi(t)=0,x<\tilde{A}(t)\le x+{\rm d}x\},\quad i=0,1,2\ldots,\\ & Q_{i,j}(x,t)\,{\rm d}x=\Pr\{N(t)=i,\xi(t)=j,x<\tilde{A}(t)\le x+{\rm d}x\},\quad i=0,1,2\ldots,\ j=1,\ldots,h. \end{align*}

Employing the above-stated probabilities, and considering the state transitions between time $t$ and $t+\Delta t$ like the usual arguments as in the birth and death model, we can have the following partial differential equations (1) to (6) for the tri-variate Markov process.

(1)\begin{align} \left(\frac{\partial}{\partial t}-\frac{\partial}{\partial x}\right)P_0(x,t)& =\sum_{j=1}^h\theta_jQ_{0,j}(x,t), \end{align}

(2)\begin{align} \left(\frac{\partial}{\partial t}-\frac{\partial}{\partial x}\right)P_n(x,t)& ={-}\mu_0P_n(x,t)+\mu_0 P_{n+1}(x,t)+a(x)\sum_{k=1}^n g_kP_{n-k}(0,t)\nonumber\\ & \quad +\sum_{j=1}^h\theta_jQ_{n,j}(x,t),\quad n=1,2,\ldots,b-1, \end{align}

(3)\begin{align} \left(\frac{\partial}{\partial t}-\frac{\partial}{\partial x}\right)P_n(x,t)& ={-}\mu_0P_n(x,t)+\mu_0 P_{n+1}(x,t)+a(x)\sum_{k=1}^bg_kP_{n-k}(0,t)\nonumber\\ & \quad +\sum_{j=1}^h\theta_jQ_{n,j}(x,t), \quad n\ge b, \end{align}

(4)\begin{align} \left(\frac{\partial}{\partial t}-\frac{\partial}{\partial x}\right)Q_{0,j}(x,t)& ={-}\theta_jQ_{0,j}(x,t)+\mu_jQ_{1,j}(x,t)+\mu_0\alpha _jP_1(x,t),\quad j=1,2,\ldots,h, \end{align}

(5)\begin{align} \left(\frac{\partial}{\partial t}-\frac{\partial}{\partial x}\right)Q_{n,j}(x,t)& ={-}(\theta_j+\mu_j)Q_{n,j}(x,t)+\mu_jQ_{n+1,j}(x,t)\nonumber\\ & \quad +a(x)\sum_{k=1}^ng_kQ_{n-k,j}(0,t), \quad n=1,\ldots,b-1,\ j=1,\ldots,h, \end{align}

(6)\begin{align} \left(\frac{\partial}{\partial t}-\frac{\partial}{\partial x}\right)Q_{n,j}(x,t)& ={-}(\theta_j+\mu_j)Q_{n,j}(x,t)+\mu_jQ_{n+1,j}(x,t)\nonumber\\ & \quad +a(x)\sum_{k=1}^bg_kQ_{n-k,j}(0,t), \quad n\ge b,\ j=1,\ldots,h. \end{align}

As many articles have pointed out, some probabilistic interpretations will be helpful to understand how these equations are derived. For example, Eq. (6) is obtained at time $t+\Delta t$ considering all possibilities. Note that when time $t$ is increased by $\Delta t$, the remaining inter-batch arrival time will be reduced by $x-\Delta t$. Thus, we have

\begin{align*} Q_{n,j}(x-\Delta t, t+\Delta t)& =Q_{n,j}(x,t)(1-\theta_j\Delta t+o(\Delta t))(1-\mu_j\Delta t+o(\Delta t))+\mu_j\Delta tQ_{n+1,j}(x,t)\\ & \quad +a(x)\Delta t\sum_{k=1}^bg_kQ_{n-k,j}(0,t),\quad n\ge b,\ j=1,\ldots,h. \end{align*}

The above equation explains possible cases for the probability that there are $n(n\ge b)$ customers in the system and the server is on type $j$ vacation when the remaining inter-batch arrival time is $x-\Delta t$ at time $t+\Delta t$. In a similar way, Eqs. (1)–(5) can be also obtained.

Let $\lim _{t\rightarrow \infty }P_n(x,t)=P_n(x)$ and $\lim _{t\rightarrow \infty }Q_{n,j}(x,t)=Q_{n,j}(x)$, the Kolmogorov forward equations governing the system in steady-state for the proposed model is:

(7)\begin{align} -\frac{\rm d}{{\rm d}x}P_0(x)& =\sum_{j=1}^h\theta_jQ_{0,j}(x), \end{align}

(8)\begin{align} -\frac{\rm d}{{\rm d}x}P_n(x)& ={-}\mu_0P_n(x)+\mu_0 P_{n+1}(x)+a(x)\sum_{k=1}^ng_kP_{n-k}(0)\nonumber\\ & \quad +\sum_{j=1}^h \theta_jQ_{n,j}(x),\quad n=1,\ldots,b-1, \end{align}

(9)\begin{align} -\frac{\rm d}{{\rm d}x}P_n(x)& ={-}\mu_0P_n(x)+\mu_0P_{n+1}(x)+a(x)\sum_{k=1}^bg_kP_{n-k}(0)\nonumber\\ & \quad+\sum_{j=1}^h \theta_jQ_{n,j}(x),\quad n\ge b, \end{align}

(10)\begin{align} -\frac{\rm d}{{\rm d}x}Q_{0,j}(x)& ={-}\theta_jQ_{0,j}(x)+\mu_jQ_{1,j}(x)+\mu_0\alpha_jP_1(x),\quad j=1,2,\ldots,h, \end{align}

(11)\begin{align} -\frac{\rm d}{{\rm d}x}Q_{n,j}(x)& ={-}(\theta_j+\mu_j)Q_{n,j}(x)+\mu_jQ_{n+1,j}(x)\nonumber\\ & \quad +a(x)\sum_{k=1}^n g_kQ_{n-k,j}(0),\quad n=1,\ldots,b-1,\ j=1,\ldots,h, \end{align}

(12)\begin{align} -\frac{\rm d}{{\rm d}x}Q_{n,j}(x)& ={-}(\theta_j+\mu_j)Q_{n,j}(x)+\mu_jQ_{n+1,j}(x)\nonumber\\ & \quad+a(x)\sum_{k=1}^b\!g_kQ_{n-k,j}(0),\quad n\ge b,\ j=1,\ldots,h. \end{align}

Since the Laplace transform turns the operation of differentiation into the algebraic operation in $s$-domain, and can greatly simplify the solution of problems involving differential equations, we further define the Laplace transforms of $P_n(x)$ and $Q_{n,j}(x)$ as follows

$$P_n^*(s)=\int_0^\infty e^{{-}sx}P_n(x)\,{\rm d}x,\quad n\ge 0, \quad Q_{n,j}^*(s)=\int_0^\infty e^{{-}sx}Q_{n,j}(x)\,{\rm d}x,\quad n\ge 0,\ j=1,\ldots,h.$$

Additionally, for $s=0$, we define

$$P_n^*(0)\equiv P_n=\int_0^\infty P_n(x)\,{\rm d}x,\quad n\ge 0,\quad Q_{n,j}^*(0)\equiv Q_{n,j}=\int_0^\infty Q_{n,j}(x)\,{\rm d}x,\quad n\ge 0,\ j=1,\ldots,h.$$

Thus, $P_n$ is the probability that there are $n$ customers in the system when the server is in a normal busy period. Similarly, $Q_{n,j}$ is the probability that there are $n$ customers in the system while the server is on type $j$ vacation. Taking the Laplace transform on both sides of Eqs. (7)–(12), the corresponding algebraic equations are given by

(13)\begin{align} &sP_0^*(s)=P_0(0)-\sum_{j=1}^h\theta_jQ_{0,j}^*(s), \end{align}

(14)\begin{align} & (s-\mu_0)P_n^*(s)+\mu_0P_{n+1}^*(s)=P_n(0)-a^*(s)\sum_{k=1}^n g_kP_{n-k}(0)\nonumber\\ & \quad -\sum_{j=1}^h\theta_jQ_{n,j}^*(s),\quad n=1,\ldots,b-1, \end{align}

(15)\begin{align} & (s-\mu_0)P_n^*(s)+\mu_0P_{n+1}^*(s)=P_n(0)-a^*(s)\sum_{k=1}^b g_kP_{n-k}(0)-\sum_{j=1}^h \theta_jQ_{n,j}^*(s),\quad n\ge b, \end{align}

(16)\begin{align} & (s-\theta_j)Q_{0,j}^*(s)+\mu_jQ_{1,j}^*(s)=Q_{0,j}(0)-\mu_0\alpha_jP_1^*(s),\quad j=1,\ldots,h, \end{align}

(17)\begin{align} & [s-(\theta_j+\mu_j)]Q_{n,j}^*(s)+\mu_jQ_{n+1,j}^*(s)=Q_{n,j}(0)-a^*(s)\sum_{k=1}^n g_kQ_{n-k,j}(0),\nonumber\\ & \qquad n=1,\ldots,b-1,\ j=1,\ldots,h, \end{align}

(18)\begin{align} & [s-(\theta_j+\mu_j)]Q_{n,j}^*(s)+\mu_jQ_{n+1,j}^*(s)=Q_{n,j}(0)-a^*(s)\sum_{k=1}^b g_kQ_{n-k,j}(0),\nonumber\\ & \qquad n\ge b,\quad j=1,\ldots,h. \end{align}

Adding Eqs. (13) to (18), term by term on both sides, it yields

$$\sum_{n=0}^\infty P_n^*(s)+\sum_{n=0}^\infty\sum_{j=1}^hQ_{n,j}^*(s)=\frac{1-a^*(s)}{s}\left(\sum_{n=0}^\infty P_n(0)+\sum_{n=0}^\infty\sum_{j=1}^hQ_{n,j}(0)\right).$$

Notice that $P_n^*(s)\rightarrow P_n$, $Q_{n,j}^*(s)\rightarrow Q_{n,j}$ and $1-a^*(s)\rightarrow 0$ as $s\rightarrow 0$, so L'H$\hat {{\rm o}}$spital's Rule and the normalization condition give

$$1=\sum_{n=0}^\infty P_n+\sum_{n=0}^\infty\sum_{j=1}^hQ_{n,j}=\frac{1}{\lambda}\left(\sum_{n=0}^\infty P_n(0)+\sum_{n=0}^\infty\sum_{j=1}^hQ_{n,j}(0)\right).$$

Let $P_n^-$ denote the probability that an arriving batch finds the server in a normal working mode and sees $n$ customers in the system, and $Q_{n,j}^-$ is defined as the probability that an incoming batch of the customers also finds $n$ customers in the system and the server is on type $j$ vacation. According to Bayes’ theorem, the following formulas are obtained:

(19)\begin{align} P_n^-& =\Pr\{n \ \text{customers in the system prior to an arrival of a batch when the server}\nonumber\\ & \qquad \text{is in normal working mode}\,|\,\text{a group of customers will arrive soon}\}\nonumber\\ & =\frac{P_n(0)}{\sum_{n=0}^\infty P_n(0)+\sum_{n=0}^\infty\sum_{j=1}^hQ_{n,j}(0)}=\frac{1}{\lambda}P_n(0),\quad n=0,1,\ldots, \end{align}

(20)\begin{align} Q_{n,j}^-& =\Pr\{n\ \text{customers in the system prior to an arrival of a batch when the server}\nonumber\\ & \qquad \text{is on type}\ j\ \text{vacation}\,|\,\text{a group of customers will arrive soon}\}\nonumber\\ & =\frac{Q_{n,j}(0)}{\sum_{n=0}^\infty P_n(0)+\sum_{n=0}^\infty\sum_{j=1}^hQ_{n,j}(0)}=\frac{1}{\lambda}Q_{n,j}(0),\quad n=0,1,\ldots,\ j=1,\ldots,h. \end{align}

Clearly, it may be seen that once the expressions of $P_n(0)$, $P_n^*(s)$, $Q_{n,j}(0)$, and $Q_{n,j}^*(s)$ are given, the queue-length distribution both at pre-arrival ($P_n^-$ and $Q_{n,j}^-$) and arbitrary epochs ($P_n$ and $Q_{n,j}$) can be determined from these quantities. We will address this topic in the following section by using the shift operator method for solving the sets of difference equations (see (13)–(18)) that arise in analyzing such a queue.

4. Stationary queue-length distributions at two different epochs

For the purpose of analysis, the discrete variable $n$ which takes on values in nonnegative integers will be viewed as the independent variable while $P_n(0)$, $P_n^*(s)$, $Q_{n,j}(0)$, and $Q_{n,j}^*(s)$ (whose value “depends” on the value of $n$) can be viewed as functions of $n$. The forward shift operator $\mathcal {D}$ acting on the sequences $\{P_n(0),n\ge 0\}$, $\{P_n^*(s),n\ge 0\}$, $\{Q_{n,j}(0), n\ge 0\}$, and $\{Q_{n,j}^*(s), n\ge 0\}$ is defined by

\begin{align*} & \mathcal{D}^lP_n(0)=P_{n+l}(0),\quad \mathcal{D}^lP_n^*(s)=P_{n+l}^*(s),\quad l\ge 1,\\ & \mathcal{D}^lQ_{n,j}(0)=Q_{n+l,j}(0),\quad \mathcal{D}^lQ_{n,j}^*(s)=Q_{n+l,j}^*(s),\quad l\ge 1,\ j=1,\ldots,h. \end{align*}

With the aid of the notation of the forward shift operator $\mathcal {D}$, the difference equation (18) can be written as

(21)\begin{equation} [s-(\theta_j+\mu_j)+\mu_j\mathcal{D}]Q_{n,j}^*(s)= \left(\mathcal{D}^b-a^*(s)\sum_{k=1}^bg_k\mathcal{D}^{b-k}\right)Q_{n-b,j}(0), \quad n\ge b,\ j=1,\ldots,h. \end{equation}

Re-indexing the variable of $Q_{n-b,j}(0)$ as $n-b\rightarrow n$, and setting $s=\theta _j+\mu _j-\mu _j\mathcal {D}$ for $j=1,\ldots,h$, Eq. (21) can be reduced to a homogeneous difference equation with constant coefficients

(22)\begin{equation} \left(\mathcal{D}^b-a^*(\theta_j+\mu_j-\mu_j\mathcal{D})\sum_{k=1}^bg_k\mathcal{D}^{b-k}\right)Q_{n,j}(0)=0,\quad n\ge 0,\ j=1,\ldots,h. \end{equation}

According to the basic theory of difference equation, for a fixed $j$, the characteristic equation associated with Eq. (22) is

(23)\begin{equation} z^b-a^*(\theta_j+\mu_j-\mu_jz)\sum_{k=1}^bg_kz^{b-k}=0,\quad j=1,\ldots,h. \end{equation}

Next, we use Rouch$\acute {\rm e}$'s theorem to find the number of zeros of the function $p(z)=z^b-a^*(\theta _j+\mu _j-\mu _jz)\sum _{k=1}^bg_kz^{b-k}$ that lies inside the unit circle. Our results will be presented in the form of the following lemma.

Lemma 1. If $\lambda \overline {g}/\mu _j<1$, then $p(z)$ has $b$ zeros (counted with multiplicity) in the disk $\{|z|<1\}$.

Proof. To apply Rouch$\acute {\rm e}$'s theorem, we seek to express $p(z)$ in the form $p(z)=f(z)+h(z)$, where the function $f(z)$ dominates $h(z)$ on the unit circle, and where it is apparent how many zeros has inside the unit circle. Thus, our choice for $f(z)$ in this case is $f(z)=z^b$, which has $b$ zeros inside the unit circle, all at the origin. Furthermore, we take $h(z)=-a^*(\theta _j+\mu _j-\mu _jz)\sum _{k=1}^bg_kz^{b-k}$. For notational convenience, let $U(z)=a^*(\theta _j+\mu _j-\mu _jz)$. When $\epsilon >0$ is sufficiently small, $U(z)$ is analytic in $|z|\le 1+\epsilon$. Thus, according to the Taylor's theorem for analytic function, $U(z)$ can be represented as a power series $\sum _{l=0}^\infty u_l(z-1)^l$. Also, we have formulas for the coefficients $u_l={U^{(l)}(1)}/{l!}=({1}/{2\pi {\mathbf {i}}})\oint _\Gamma ({U(z)}/{(z-1)^{l+1}})\,{\rm d}z$, where $\mathbf {i}$ is a square root of $-1$, and $\Gamma$ is any simple closed curve in $|z|\le 1+\epsilon$ around $1$. Now, employing the Taylor expansion for $U(z)$, we may estimate $|h(z)|$ and $|f(z)|$ on the simple closed curve $|z|=1-\delta$, where $\delta$ is a small positive real number. Since

\begin{align*} |h(z)|_{z=1-\delta}& = |U(z)|\left|\sum_{k=1}^bg_kz^{b-k}\right|_{z=1-\delta}\le U(|z|)\sum_{k=1}^bg_k|z|^{b-k}_ {z=1-\delta}=U(1-\delta)\sum_{k=1}^bg_k(1-\delta)^{b-k}\\ & =\left[U(1)+\frac{U^\prime(1)}{1!}(1-\delta-1)+\sum_{l=2}^\infty\frac{U^{(l)}(1)}{l!}(1-\delta-1)^l\right] \sum_{k=1}^b g_k[1-(b-k)\delta+{\rm o}(\delta)]\\ & = \left[\int_0^\infty e^{-\theta_jt}\,{\rm d}A(t)-\delta\!\int_0^\infty e^{-\theta_jt}\mu_jt\,{\rm d}A(t)+{\rm o}(\delta)\right]\sum_{k=1}^bg_k[1-(b-k)\delta+{\rm o}(\delta)]\\ & = \left[\int_0^\infty e^{-\theta_jt}(1-\delta\mu_jt)\,{\rm d}A(t)+{\rm o}(\delta)\right]\sum_{k=1}^bg_k[1-(b-k)\delta+{\rm o}(\delta)]\\ & \le \left[\int_0^\infty(1-\delta\mu_jt)\,{\rm d}A(t)+{\rm o}(\delta)\right]\sum_{k=1}^bg_k[1-(b-k)\delta+{\rm o}(\delta)]\\ & = \left[1-\delta\frac{\mu_j}{\lambda}+{\rm o}(\delta)\right][1-b\delta+\overline{g}\delta+{\rm o}(\delta)]=1-b\delta-\left( \frac{\mu_j}{\lambda}-\overline{g}\right)\delta+{\rm o}(\delta)\\ & <|f(z)|_{z=1-\delta}=(1-\delta)^b=1-b\delta+{\rm o}(\delta), \end{align*}

and both $f(z)$ and $h(z)$ are analytic for $|z|\le 1-\delta$, letting $\delta$ tend to zero, Rouch$\acute {\rm e}$'s theorem tells us that $f(z)$ and $p(z)=f(z)+h(z)$ have the same number of zeros in $|z|<1$, which is $b$.

On the other hand, in the existing literature, authors like Chaudhry [8] and Tijms [25] have emphasized that root-finding in queueing theory is well structured, in the sense that the roots of the characteristic equation are distinct for most queueing models. Therefore, based on the above point of view, we assume that the $b$ roots of the characteristic equation (23) are distinct, denoted by $r_{m,j}$, $m=1,\ldots,b$, $j=1,\ldots,h$. Then, an immediate consequence for this particular case is that the general solution of the homogeneous equation (22) is as follows

(24)\begin{equation} Q_{n,j}(0)=\sum_{m=1}^bc_{m,j}r_{m,j}^n,\quad n\ge 0,\ j=1,\ldots,h, \end{equation}

where for fixed $j$, $c_{1,j},c_{2,j},\ldots$, and $c_{b,j}$ are real or complex constants whose values can be determined by a system of linear equations (see subsequent discussions). Substitution of Eq. (24) into Eq. (21) gives

(25)\begin{align} & [s-(\theta_j+\mu_j)+\mu_j\mathcal{D}]Q_{n,j}^*(s)=\sum_{m=1}^bc_{m,j}r_{m,j}^n-a^*(s)\sum_{k=1}^bg_k \left(\sum_{m=1}^bc_{m,j}r_{m,j}^{n-k}\right)\nonumber\\ & \quad =\sum_{m=1}^bc_{m,j}r_{m,j}^n-a^*(s)\sum_{m=1}^bc_{m,j}\sum_{k=1}^bg_kr_{m,j}^{n-k},\quad n\ge b,\ j=1,\ldots,h. \end{align}

Any solution $Q_{n,j}^*(s)$ of Eq. (25) may be written as

$$Q_{n,j}^*(s)=Q_{n,j}^{*{\rm (par)}}(s)+Q_{n,j}^{*{\rm (hom)}}(s)=Q_{n,j}^{*{\rm (par)}}(s)+\mathbb{C}_1\left(1+\frac{\theta_j-s}{\mu_j}\right)^n,\quad n\ge b,$$

where $Q_{n,j}^{*{\rm (par)}}(s)$ is a particular solution of the nonhomogeneous equation (25), and $Q_{n,j}^{*{\rm (hom)}}(s)=\mathbb {C}_1(1+{(\theta _j-s)}/{\mu _j})^n$ is the general solution of the corresponding homogenous equation $[s-(\theta _j+\mu _j)+ \mu _j\mathcal {D}]Q_{n,j}^*(s)=0$. It is easy to verify that the undetermined constant $\mathbb {C}_1=0$ because of the condition $\sum _{n=b}^\infty Q_{n,j}^*(0)=\sum _{n=b}^\infty Q_{n,j}<\infty$. Thus, for $n\ge b$, we finally have $Q_{n,j}^*(s)=Q_{n,j}^{*{\rm (par)}}(s)$.

Now, we focus our attention on finding a particular solution of Eq. (25). Since the nonhomogeneous term is a finite linear combination of the function $r_{m,j}^n$, the method of undetermined coefficient is used to find $Q_{n,j}^{*{\rm (par)}}(s)$. Here, we can guess the trial solutions with undetermined coefficients, plug them into the difference equation (25), and then solve for the unknown coefficients to obtain the particular solution as below

(26)\begin{equation} Q_{n,j}^*(s)=Q_{n,j}^{*{\rm (par)}}(s)=\sum_{m=1}^bc_{m,j}\frac{r_{m,j}^n-a^*(s)\sum_{k=1}^bg_kr_{m,j}^{n-k}}{s-\theta_j-\mu_j(1-r_{m,j})},\quad n\ge b,\ j=1,\ldots,h. \end{equation}

Since Eqs. (17) and (18) have almost exactly the same form, we now seek the suitable conditions under which $Q_{n,j}^*(s)|_{s=0}$ has the same expression as in Eq. (26) for $1\le n\le b-1$, when $s=0$. This point can be justified by the matrix geometric algorithm (see Neuts [21]), and is not our intuitive guess. Actually, by assuming that the size of an arriving batch is bounded, the jumps to the right of the embedded Markov chain are bounded (although our approach can avoid discussion of the embedded Markov chain). This implies that the Markov chain fits into the standard matrix-geometric framework by appropriately defining the levels, and its stationary distribution is readily available in matrix-geometric form. In other words, for $n=0,1,2,\ldots$, if we let $\bf \pi _n=(P^-_n,Q^-_{n,1},Q^-_{n,2},\ldots, Q^-_{n,h})$, Neuts has shown that there exists a positive matrix $\bf R$ such that $\bf \pi _n=\bf \pi _0\bf R^n$. Furthermore, from the results regarding the steady-state queue-length distributions immediately before batch arrival, we can derive the stationary queue-length distribution at arbitrary epochs by employing the classical arguments based on renewal theory and semi-Markov process. For the above reasons, and noting that $P_n^*(s)|_{s=0}=P_n$ and $Q_{n,j}^*(s)|_{s=0}=Q_{n,j}$, we think the expressions of $Q_{n,j}^*(0)$ and $P_n^*(0)$ can necessarily be unified regardless of $n\ge b$ or not. Hence, by putting Eq. (24) into Eqs. (18) and (17), and comparing the last term of the right-hand side of Eqs. (18) and (17) gives the following relation, which the $c_{m,j}$ must satisfy:

(27)\begin{align} & \sum_{k=1}^ng_k\sum_{m=1}^bc_{m,j}r_{m,j}^{n-k}=\sum_{k=1}^bg_k\sum_{m=1}^bc_{m,j}r_{m,j}^{n-k} \Rightarrow\sum_{m=1}^bc_{m,j}\sum_{k=1}^ng_kr_{m,j}^{n-k}=\sum_{m=1}^bc_{m,j} \sum_{k=1}^bg_kr_{m,j}^{n-k}\nonumber\\ & \quad \Rightarrow\sum_{m=1}^bc_{m,j}\sum_{k=n+1}^bg_kr_{m,j}^{n-k}=0,\quad n=b-1,b-2,\ldots,1,\ j=1,2,\ldots,h. \end{align}

Notice that for $g_b\neq 0$, setting $n=b-1$, $b-2$, $\ldots$, $1$ in Eq. (27), respectively, allows us to conclude that

(28)\begin{align} & \left\{\begin{array}{l} \displaystyle\sum_{m=1}^b\dfrac{c_{m,1}}{r_{m,1}}=0\\ \displaystyle\sum_{m=1}^b\dfrac{c_{m,1}}{r_{m,1}^2}=0\\ \vdots \\ \displaystyle\sum_{m=1}^b\dfrac{c_{m,1}}{r_{m,1}^{b-2}}=0\\ \displaystyle\sum_{m=1}^b\dfrac{c_{m,1}}{r_{m,1}^{b-1}}=0 \end{array}\right.,\quad \left\{\begin{array}{l} \displaystyle\sum_{m=1}^b\dfrac{c_{m,2}}{r_{m,2}}=0\\ \displaystyle\sum_{m=1}^b\dfrac{c_{m,2}}{r_{m,2}^2}=0\\ \vdots \\ \displaystyle\sum_{m=1}^b\dfrac{c_{m,2}}{r_{m,2}^{b-2}}=0\\ \displaystyle\sum_{m=1}^b\dfrac{c_{m,2}}{r_{m,2}^{b-1}}=0 \end{array}\right.,\ldots\ldots, \left\{\begin{array}{l} \displaystyle\sum_{m=1}^b\dfrac{c_{m,h-1}}{r_{m,h-1}}=0\\ \displaystyle\sum_{m=1}^b\dfrac{c_{m,h-1}}{r_{m,h-1}^2}=0\\ \vdots \\ \displaystyle\sum_{m=1}^b\dfrac{c_{m,h-1}}{r_{m,h-1}^{b-2}}=0\\ \displaystyle\sum_{m=1}^b\dfrac{c_{m,h-1}}{r_{m,h-1}^{b-1}}=0 \end{array}\right.,\quad \left\{\begin{array}{l} \displaystyle\sum_{m=1}^b\dfrac{c_{m,h}}{r_{m,h}}=0\\ \displaystyle\sum_{m=1}^b\dfrac{c_{m,h}}{r_{m,h}^2}=0\\ \vdots \\ \displaystyle\sum_{m=1}^b\dfrac{c_{m,h}}{r_{m,h}^{b-2}}=0\\ \displaystyle\sum_{m=1}^b\dfrac{c_{m,h}}{r_{m,h}^{b-1}}=0 \end{array}\right.. \end{align}

For $j=1,2,\ldots,h$, let $\bf c_j=(c_{1,j},c_{2,j},\ldots,c_{b,j})^{\mathsf {T}}$ and

$$\bf R_j=\left(\begin{array}{ccccc} r_{1,j}^{{-}1} & r_{2,j}^{{-}1} & \cdots & r_{b-1,j}^{{-}1} & r_{b,j}^{{-}1} \\ r_{1,j}^{{-}2} & r_{2,j}^{{-}2} & \cdots & r_{b-1,j}^{{-}2} & r_{b,j}^{{-}2} \\ \vdots & \vdots & \vdots & \vdots & \vdots\\ r_{1,j}^{-(b-2)} & r_{2,j}^{-(b-2)} & \cdots & r_{b-1,j}^{-(b-2)} & r_{b,j}^{-(b-2)}\\ r_{1,j}^{-(b-1)} & r_{2,j}^{-(b-1)} & \cdots & r_{b-1,j}^{-(b-1)} & r_{b,j}^{-(b-1)} \end{array}\right),$$

where “$\mathsf {T}$” represents the transpose operation on a vector or a matrix. With these notations, Eq. (27) may be written in matrix form as

(29)\begin{equation} \left(\begin{array}{cccccc} \bf R_1 & \bf 0 & \bf 0 & \cdots & \bf 0 & \bf 0 \\ \bf 0 & \bf R_2 & \bf 0 & \cdots & \bf 0 & \bf 0 \\ \bf 0 & \bf 0 & \bf R_3 & \cdots & \bf 0 & \bf 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \bf 0 & \bf 0 & \bf 0 & \cdots & \bf R_{h-1} & \bf 0 \\ \bf 0 & \bf 0 & \bf 0 & \cdots & \bf 0 & \bf R_h\\ \end{array}\right)\left(\begin{array}{c} \bf c_1\\ \bf c_2\\ \bf c_3\\ \vdots \\ \bf c_{h-1}\\ \bf c_h \end{array}\right)=\left(\begin{array}{c} \bf 0\\ \bf 0\\ \bf 0\\ \vdots \\ \bf 0\\ \bf 0 \end{array}\right), \end{equation}

where $\bf 0$ denotes a zero matrix or a zero column vector of appropriate dimension. If the constants $c_{1,1}$, $\ldots$, $c_{b,1}$, $\ldots$, $c_{1,h}$, $\ldots$, $c_{b,h}$ satisfy Eq. (28), the same expression for $Q_{n,j}^*(s)$ can be derived from Eqs. (17) and (18). This also means that for any $n\ge 1$, we have

(30)\begin{equation} Q_{n,j}^*(s)=\sum_{m=1}^bc_{m,j}\frac{r_{m,j}^n-a^*(s)\sum_{k=1}^bg_kr_{m,j}^{n-k}}{s-\theta_j-\mu_j(1-r_{m,j})},\quad j=1,\ldots,h. \end{equation}

However, since the number of unknowns is greater than the number of equations, the constants $c_{m,j}$, $m=1,\ldots,b$, $j=1,\ldots,h$ can not be uniquely determined by Eq. (29). So we need to find more conditions to get $c_{m,j}$.

Again, as in Eq. (18) treated above, Eq. (15) can also be expressed in terms of the forward shift operator $\mathcal {D}$ as below

(31)\begin{equation} (s-\mu_0+\mu_0\mathcal{D})P_n^*(s)=\left(\mathcal{D}^b-a^*(s)\sum_{k=1}^b\!g_k\mathcal{D}^{b-k}\right)P_{n-b}(0)-\sum_{j=1}^h \theta_jQ_{n,j}^*(s),\quad n\ge b. \end{equation}

Letting $s=\mu _0-\mu _0\mathcal {D}$ on both sides of Eq. (30) gives

(32)\begin{equation} \left(\mathcal{D}^b-a^*(\mu_0-\mu_0\mathcal{D})\sum_{k=1}^b\!g_k\mathcal{D}^{b-k}\right)P_{n-b}(0) =\sum_{j=1}^h\theta_jQ_{n,j}^*(\mu_0-\mu_0\mathcal{D}),\quad n\ge b. \end{equation}

Similarly, by setting $n\ge 0$ instead of $n\ge b$ and substituting Eq. (30) into Eq. (32) yields

(33)\begin{align} & \left(\mathcal{D}^b-a^*(\mu_0-\mu_0\mathcal{D})\sum_{k=1}^bg_k\mathcal{D}^{b-k}\right)P_n(0) =\sum_{j=1}^h\theta_jQ_{n+b,j}^*(\mu_0-\mu_0\mathcal{D})\nonumber\\ & \quad =\sum_{j=1}^h\theta_j\sum_{m=1}^bc_{m,j}\left\{\frac{r_{m,j}^b-a^*(\mu_0-\mu_0\mathcal{D})\sum_{k=1}^bg_kr_{m,j}^{b-k}} {\mu_0-\mu_0\mathcal{D}-\theta_j-\mu_j(1-r_{m,j})}\right\}r_{m,j}^n,\quad n\ge 0. \end{align}

Eq. (33) is a nonhomogeneous equation for $P_n(0)$ as a function of $n$. The general solution of Eq. (33) consists of the sum of the solution to the homogeneous equation

(34)\begin{equation} \left(\mathcal{D}^b-a^*(\mu_0-\mu_0\mathcal{D})\sum_{k=1}^bg_k\mathcal{D}^{b-k}\right)P_n(0)=0,\quad n\ge 0, \end{equation}

and any particular solution of Eq. (33). Notice that the characteristic equation associated with Eq. (34) is

(35)\begin{equation} z^b-a^*(\mu_0-\mu_0z)\sum_{k=1}^bg_kz^{b-k}=0. \end{equation}

If we start with Eq. (35) and still choose $f(z)=z^b$ and $q(z)=-a^*(\mu _0-\mu _0z)\sum _{k=1}^b g_kz^{b-k}$, analogous to the proof of Lemma 1, we can prove the characteristic equation (35) has exactly $b$ roots inside the unit disk under the assumption that $\overline {g}\lambda <\mu _0$. Let these roots be denoted by $\omega _1$, $\omega _2$, $\ldots$, $\omega _b$, and also assume that they are distinct. Then the general solution to the homogeneous equation (34) is $\sum _{m=1}^bf_m\omega _m^n$, where $f_1$, $f_2$, $\ldots$, and $f_b$ are $b$ arbitrary constants that need to be determined. Next, to find a particular solution of Eq. (33), we observe the structure of the nonhomogeneous term, and guess that it has to have the form

(36)\begin{equation} P_n^{({\rm par})}(0)=\sum_{j=1}^h\sum_{m=1}^bd_{m,j}r_{m,j}^n,\quad n\ge0, \end{equation}

where $d_{m,j}$ can be expressed in terms of $c_{m,j}$ and $r_{m,j}$ by substituting Eq. (36) into Eq. (33). Doing this gives

$$d_{m,j}=\frac{c_{m,j}\theta_j}{(\mu_0-\mu_j)(1-r_{m,j})-\theta_j},\quad m=1,\ldots,b,\ j=1,\ldots,h.$$

Hence final general solution to the nonhomogeneous equation (32) is given by

(37)\begin{equation} P_n(0)=\sum_{m=1}^bf_m\omega_m^n+\sum_{j=1}^h\sum_{m=1}^b \frac{c_{m,j}\theta_j}{(\mu_0-\mu_j)(1-r_{m,j})-\theta_j}r_{m,j}^n, \quad n\ge 0. \end{equation}

The task now is to find the expression of $P_n^*(s)$. Plugging Eqs. (30) and (37) into Eq. (31), and after some algebraic manipulation, we get

(38)\begin{align} & (s-\mu_0+\mu_0\mathcal{D})P_n^*(s)\nonumber\\ & \quad =\left(\mathcal{D}^b-a^*(s)\sum_{k=1}^bg_k\mathcal{D}^{b-k}\right)\left( \sum_{m=1}^b f_m\omega_m^{n-b}+\sum_{j=1}^h\sum_{m=1}^b \frac{c_{m,j}\theta_jr_{m,j}^{n-b}}{(\mu_0-\mu_j)(1-r_{m,j}) -\theta_j}\right)\nonumber\\ & \qquad -\sum_{j=1}^h\theta_j\sum_{m=1}^bc_{m,j}\frac{r_{m,j}^n-a^*(s)\sum_{k=1}^bg_kr_{m,j}^{n-k}} {s-\theta_j-\mu_j(1-r_{m,j})}\nonumber\\ & \quad =\sum_{m=1}^bf_m\left(\omega_m^n-a^*(s)\sum_{k=1}^b g_k\omega_m^{n-k}\right)+\sum_{j=1}^h \sum_{m=1}^b \frac{c_{m,j}\theta_j(r_{m,j}^n-a^*(s) \sum_{k=1}^bg_kr_{m,j}^{n-k})}{(\mu_0-\mu_j)(1-r_{m,j})-\theta_j}\nonumber\\ & \qquad -\sum_{j=1}^h\sum_{m=1}^bc_{m,j}\theta_j\left\{\frac{r_{m,j}^n-a^*(s)\sum_{k=1}^bg_kr_{m,j}^{n-k}} {s-\theta_j-\mu_j(1-r_{m,j})}\right\}\nonumber\\ & \quad =\sum_{m=1}^b f_m\left(1-a^*(s)\sum_{k=1}^b g_k\omega_m^{{-}k}\right)\omega_m^n\nonumber\\ & \qquad +\sum_{j=1}^h\sum_{m=1}^b \frac{c_{m,j}\theta_j(1-a^*(s)\sum_{k=1}^bg_kr_{m,j}^{{-}k})[s-\mu_0(1-r_{m,j})]} {[(\mu_0-\mu_j)(1-r_{m,j})-\theta_j][s-\theta_j-\mu_j(1-r_{m,j})]}r_{m,j}^n,\quad n\ge b. \end{align}

This indicates that the sequence $\{P_n^*(s),n\ge b\}$ satisfies the first-order nonhomogeneous difference equation (38). Similar to our previous discussion, we of course first find the general solution to the corresponding homogeneous equation of Eq. (38). Here, the general solution is given as $P_n^{*({\rm hom})}(s)=\mathbb {C}_2(1-{s}/{\mu _0})^n$, and $\mathbb {C}_2$ is an undetermined constant. Notice also that the nonhomogeneous term is a linear combination of $\omega _m^n$ and $r_{m,j}^n$, a particular solution of Eq. (38) can be given by

(39)\begin{align} P_n^{*({\rm par})}(s)& =\sum_{m=1}^b\frac{f_m[1-a^*(s)G(\omega_m^{{-}1})]}{s-\mu_0(1-\omega_m)} \omega_m^n\nonumber\\ & \quad +\sum_{j =1}^h\sum_{m=1}^b\frac{c_{m,j}\theta_j[1-a^*(s)G(r_{m,j}^{{-}1})]r_{m,j}^n} {[(\mu_0-\mu_j)(1-r_{m,j})-\theta_j][s-\theta_j-\mu_j(1-r_{m,j})]}, \quad n\ge b. \end{align}

Thus, for $n\ge b$, the general solution to Eq. (38) is the sum of the homogeneous solution and the particular solution (see Eq. (39)),

(40)\begin{align} P_n^*(s)& =\mathbb{C}_2\left(1-\frac{s}{\mu}\right)^n+\sum_{m=1}^b \frac{f_m[1-a^*(s)G(\omega_m^{{-}1})]}{s-\mu_0(1-\omega_m)}\omega_m^n\nonumber\\ & \quad +\sum_{j=1}^h\sum_{m=1}^b\frac{c_{m,j}\theta_j[1-a^*(s)G(r_{m,j}^{{-}1})]}{[(\mu_0-\mu_j)(1-r_{m,j})- \theta_j][s-\theta_j-\mu_j(1-r_{m,j})]}r_{m,j}^n,\quad n\ge b. \end{align}

Since $\sum _{n=b}^\infty P_n^*(0)=\sum _{n=b}^\infty P_n<\infty$, we must have $\mathbb {C}_2=0$. Hence, Eq. (40) reduces to

(41)\begin{align} P_n^*(s)& =\sum_{m=1}^b\frac{f_m[1-a^*(s)G(\omega_m^{{-}1}) ]}{s-\mu_0(1-\omega_m)}\omega_m^n\nonumber\\ & \quad +\sum_{j=1}^h\sum_{m=1}^b \frac{c_{m,j}\theta_j[1-a^*(s)G(r_{m,j}^{{-}1})]r_{m,j}^n}{[(\mu_0-\mu_j)(1-r_{m,j})- \theta_j][s-\theta_j-\mu_j(1-r_{m,j})]},\quad n\ge b. \end{align}

Then, we do wish to find the right conditions such that the expression for $P_n^*(s)|_{s=0}$ presented in Eq. (41) still holds when $1\le n\le b-1$ and $s=0$. We compare the second term on the right-hand side of Eqs. (14) and (15) and conclude that $f_m$ and $c_{m,j}$ satisfy the relation

(42)\begin{align} & \sum_{m=1}^b f_m \sum_{k=n+1}^b g_k\omega_m^{n-k}+\sum_{j=1}^h\sum_{m=1}^b c_{m,j} \sum_{k=n+1}^b g_k\nonumber\\ & \quad \times \frac{\theta_jr_{m,j}^{n-k}}{(\mu_0-\mu_j)(1-r_{m,j})-\theta_j}=0, \quad n=b-1,b-2,\ldots,1. \end{align}

Let us take $n=b-1,b-2,\ldots$, and $1$ in Eq. (42), respectively, and note that for $g_b\neq 0$, a system of linear equations in $b+hb$ variables can be derived from Eq. (42)

(43)\begin{equation} \left\{\begin{array}{l} \displaystyle\sum_{m=1}^b\dfrac{f_m}{\omega_m}+\sum_{j=1}^h\sum_{m=1}^b \dfrac{c_{m,j}\theta_j}{[(\mu_0-\mu_j)(1-r_{m,j})-\theta_j]r_{m,j}}=0\\ \displaystyle\sum_{m=1}^b\dfrac{f_m}{\omega_m^2}+\sum_{j=1}^h\sum_{m=1}^b \dfrac{c_{m,j}\theta_j}{[(\mu_0-\mu_j)(1-r_{m,j})-\theta_j]r_{m,j}^2}=0\\ \cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\\ \displaystyle\sum_{m=1}^b\dfrac{f_m}{\omega_m^{b-1}}+\sum_{j=1}^h\sum_{m=1}^b \dfrac{c_{m,j}\theta_j}{[(\mu_0-\mu_j)(1-r_{m,j})-\theta_j]r_{m,j}^{b-1}}=0 \end{array}\right.. \end{equation}

In other words, when $1\le n\le b-1$, if the above Eq. (43) holds, then $P_n^*(s)$ has the following uniform expression

(44)\begin{align} P_n^*(s)& = \sum_{m=1}^b\frac{f_m[1-a^*(s)G(\omega_m^{{-}1}) ]}{s-\mu_0(1-\omega_m)}\omega_m^n\nonumber\\ & \quad +\sum_{j=1}^h\sum_{m=1}^b \frac{c_{m,j}\theta_j[1-a^*(s)G(r_{m,j}^{{-}1})]r_{m,j}^n} {[(\mu_0-\mu_j)(1-r_{m,j})-\theta_j][s-\theta_j-\mu_j(1-r_{m,j})]}, \quad n\ge 1. \end{align}

Furthermore, letting $s=\theta _j$ in Eq. (16) and substituting Eqs. (24), (30), and (44) into Eq. (16) gives another linear system with $h$ equations and $hb+b$ variables.

(45)\begin{equation} \left\{\begin{array}{l} \displaystyle\sum_{m=1}^bc_{m,1}\left(\dfrac{1-a^*(\theta_1)G(r_{m,1}^{{-}1})r_{m,1}}{r_{m,1}-1}+\dfrac{\alpha_1\mu_0\theta_1[1-a^*(\theta_1)G(r_{m,1}^{- 1})]r_{m,1}}{[(\mu_0-\mu_1)(1-r_{m,1})-\theta_1]\mu_1(r_{m,1}-1)}\right)\\ \displaystyle\quad +\sum_{l=2}^h\sum_{m=1}^b\dfrac{\alpha_1\mu_0c_{m,l}\theta_l[1-a^*(\theta_1)G(r_{m,l}^{{-}1})]r_{m,l}}{[(\mu_0-\mu_l)(1- r_{m,l})-\theta_l][\theta_1-\theta_l-\mu_l(1-r_{m,l})]}\\ \displaystyle\quad +\sum_{m=1}^b\dfrac{\alpha_1\mu_0f_m[1-a^*(\theta_1) G(\omega_m^{{-}1})]\omega_m}{\theta_1-\mu_0(1-\omega_m)}=0\\ \displaystyle\sum_{m=1}^bc_{m,2}\left(\dfrac{1-a^*(\theta_2)G(r_{m,2}^{{-}1})r_{m,2}}{r_{m,2}-1}+\dfrac{\alpha_2\mu_0\theta_2[1-a^*(\theta_2)G(r_{m,2}^{- 1})]r_{m,2}}{[(\mu_0-\mu_2)(1-r_{m,2})-\theta_2]\mu_2(r_{m,2}-1)}\right)\\ \displaystyle\quad +\sum_{\substack{l=1\\ l\ne 2}}^h\sum_{m=1}^b\dfrac{\alpha_2\mu_0c_{m,l}\theta_l[1-a^*(\theta_2)G(r_{m,l}^{{-}1})]r_{m,l}}{[(\mu_0-\mu_l)(1- r_{m,l})-\theta_l][\theta_2-\theta_l-\mu_l(1-r_{m,l})]}\\ \displaystyle\quad +\sum_{m=1}^b\dfrac{\alpha_2\mu_0f_m[1-a^*(\theta_2) G(\omega_m^{{-}1})]\omega_m}{\theta_2-\mu_0(1-\omega_m)}=0\\ \cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots\\ \displaystyle\sum_{m=1}^bc_{m,h}\left(\dfrac{1-a^*(\theta_h)G(r_{m,h}^{{-}1})r_{m,h}}{r_{m,h}-1} +\dfrac{\alpha_h\mu_0\theta_h[1-a^*(\theta_h) G(r_{m,h}^{{-}1})]r_{m,h}}{[(\mu_0-\mu_h)(1-r_{m,h})-\theta_h]\mu_h(r_{m,h}-1)}\right)\\ \displaystyle\quad +\sum_{l=1}^{h-1}\sum_{m=1}^b\dfrac{\alpha_h\mu_0c_{m,l}\theta_l[1-a^*(\theta_h)G(r_{m,l}^{{-}1})]r_{m,l}}{[(\mu_0-\mu_l) (1-r_{m,l})-\theta_l][\theta_h-\theta_l-\mu_l(1-r_{m,l})]}\\ \displaystyle\quad +\sum_{m=1}^b \dfrac{\alpha_h\mu_0f_m[1-a^*(\theta_h) G(\omega_m^{{-}1})]\omega_m}{\theta_h-\mu_0(1-\omega_m)}=0\\ \end{array}\right.. \end{equation}

To enhance the calculation and simplify the programs coded in Matlab, we rewrite Eqs. (43) and (45) in the matrix form as below

(46)\begin{equation} \left(\begin{array}{ccccc} \bf H_1 & \bf H_2 & \cdots & \bf H_h & \bf W\\ \bf\Lambda_1 & \bf\Lambda_2 & \cdots & \bf\Lambda_h & \bf B \end{array}\right)\left(\begin{array}{c} \bf c_1\\ \vdots\\ \bf c_h\\ \bf f \end{array}\right)=\left(\begin{array}{c} 0\\ \vdots\\ 0\\ 0 \end{array}\right)_{(b-1+h)\times 1}, \end{equation}

where $\bf f=(f_1,f_2,\ldots,f_b)^{\mathsf {T}}$,

\begin{align*} & \bf H_j=\left(\begin{array}{cccc} \zeta_{1,j}r_{1,j}^{{-}1} & \zeta_{2,j}r_{2,j}^{{-}1} & \cdots & \zeta_{b,j}r_{b,j}^{{-}1} \\ \zeta_{1,j}r_{1,j}^{{-}2} & \zeta_{2,j}r_{2,j}^{{-}2} & \cdots & \zeta_{b,j}r_{b,j}^{{-}2} \\ \vdots & \vdots & \vdots & \vdots \\ \zeta_{1,j}r_{1,j}^{-(b-1)} & \zeta_{2,j}r_{2,j}^{-(b-1)} & \cdots & \zeta_{b,j}r_{b,j}^{-(b-1)} \end{array}\right),\quad j=1,2,\ldots,h,\\ & \bf\Lambda_j=j{\rm{th\ row}}\left(\begin{array}{ccccc} \eta_{1,j}^{(1)} & \eta_{2,j}^{(1)} & \cdots & \eta_{b-1,j}^{(1)} & \eta_{b,j}^{(1)}\\ \eta_{1,j}^{(2)} & \eta_{2,j}^{(2)} & \cdots & \eta_{b-1,j}^{(2)} & \eta_{b,j}^{(2)}\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ \psi_{1,j} & \psi_{2,j} & \cdots & \psi_{b-1,j} & \psi_{b,j}\\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \eta_{1,j}^{(h)} & \eta_{2,j}^{(h)} & \cdots & \eta_{b-1,j}^{(h)} & \eta_{b,j}^{(h)} \end{array}\right),\quad j=1,2,\ldots,h, \\ & \bf W=\left(\begin{array}{cccc} \omega_1^{{-}1} & \omega_2^{{-}1} & \cdots & \omega_b^{{-}1}\\ \omega_1^{{-}2} & \omega_2^{{-}2} & \cdots & \omega_b^{{-}2}\\ \vdots & \vdots & \vdots & \vdots \\ \omega_1^{-(b-1)} & \omega_2^{-(b-1)} & \cdots & \omega_b^{-(b-1)} \end{array}\right),\quad \bf B=\left(\begin{array}{cccc} \beta_{1,1} & \beta_{2,1} & \cdots & \beta_{b,1} \\ \beta_{1,2} & \beta_{2,2} & \cdots & \beta_{b,2} \\ \vdots & \vdots & \vdots & \vdots \\ \beta_{1,h} & \beta_{2,h} & \cdots & \beta_{b,h} \end{array}\right). \end{align*}

The elements in matrices $\bf H_j$, $\bf \Lambda _j$, and $\bf B$ are given, respectively, by

\begin{align*} \zeta_{m,j}& =\frac{\theta_j}{(\mu_0-\mu_j)(1-r_{m,j})-\theta_j},\quad m=1,\ldots,b,\ j=1,\ldots,h,\\ \eta_{m,j}^{(i)}& =\frac{\alpha_i\mu_0\theta_j[1-a^*(\theta_i)G(r_{m,j}^{{-}1})]r_{m,j}}{[(\mu_0-\mu_j) (1-r_{m,j})-\theta_j][\theta_i-\theta_j-\mu_j(1-r_{m,j})]},\quad m=1,\ldots,b, \ i,j=1,\ldots,h,\ i\neq j,\\ \psi_{m,j}& =\frac{1-a^*(\theta_j)G(r_{m,j}^{{-}1})r_{m,j}}{r_{m,j}-1}+\frac{\alpha_j\mu_0\theta_j[1-a^*(\theta_j)G(r_{m,j}^{- 1})]r_{m,j}}{[(\mu_0-\mu_j)(1-r_{m,j})-\theta_j]\mu_j(r_{m,j}-1)},\quad m=1,\ldots,b, \ j=1,\ldots,h,\\ \beta_{m,j}& =\frac{\alpha_j\mu_0[1-a^*(\theta_j) G(\omega_m^{{-}1})]\omega_m}{\theta_j-\mu_0(1-\omega_m)},\quad m=1,\ldots,b, \ j=1,\ldots,h. \end{align*}

Since $\sum _{n=0}^\infty P_n(0)+\sum _{n=0}^\infty \sum _{j=1}^hQ_{n,j}(0)=\lambda$, for $m=1,\ldots,b$ and $j=1,\ldots,h$, we can further derive a relationship for $c_{m,j}$ and $f_m$

(47)\begin{equation} \sum_{m=1}^b\frac{f_m}{1-\omega_m}+\sum_{j=1}^h\sum_{m=1}^b \frac{c_{m,j}(\mu_0-\mu_j)}{(\mu_0-\mu_j)(1-r_{m,j})- \theta_j}=\lambda. \end{equation}

Obviously, Eqs. (29), (46) and (47) constitute a system of $hb+b$ linear equations with $hb+b$ unknowns to compute the undetermined constants. Having found $c_{m,j}$ and $f_m$, we are now able to give the stationary queue-length distributions at pre-arrival and arbitrary epochs. According to Eqs. (19) and (20), we have

(48)\begin{align} & P_n^-{=}\frac{1}{\lambda}P_n(0)=\frac{1}{\lambda}\left(\sum_{m=1}^bf_m\omega_m^n +\sum_{j=1}^h\sum_{m=1}^b\frac{c_{m,j}\theta_j} {(\mu_0-\mu_j)(1-r_{m,j})-\theta_j}r_{m,j}^n\right),\quad n\ge 0, \end{align}

(49)\begin{align} & Q_{n,j}^-{=}\frac{1}{\lambda}Q_{n,j}(0)=\frac{1}{\lambda}\sum_{m=1}^bc_{m,j}r_{m,j}^n,\quad n\ge 0,\ j=1,\ldots,h. \end{align}

Taking $s=0$ in Eqs. (30) and (44), we further have

(50)\begin{align} P_n& =P_n^*(0)=\sum_{m=1}^b \frac{f_m[G(\omega_m^{{-}1})-1]\omega_m^n}{\mu_0(1-\omega_m)}\nonumber\\ & \quad +\sum_{j=1}^h\sum_{m=1}^b\frac{c_{m,j}\theta_j[G(r_{m,j}^{{-}1})-1]r_{m,j}^n}{[(\mu_0-\mu_j) (1-r_{m,j})-\theta_j][\theta_j+\mu_j(1-r_{m,j})]},\quad n\ge 1 \end{align}

(51)\begin{align} Q_{n,j}& =Q_{n,j}^*(0)=\sum_{m=1}^b\frac{c_{m,j}[G(r_{m,j}^{{-}1})-1]r_{m,j}^n}{\theta_j+\mu_j(1-r_{m,j})},\quad n\ge 1,\ j=1,\ldots,h. \end{align}

For $n=0$, $Q_{0,j}$ can be derived from Eq. (16) by making the substitution $s=0$. With some algebraic manipulation, it can be shown that

(52)\begin{align} Q_{0,j}& =Q_{0,j}^*(0)= \frac{1}{\theta_j}\left[\mu_j \sum_{m=1}^b\frac{c_{m,j}[G(r_{m,j}^{{-}1})-1] r_{m,j}}{\theta_j+\mu_j(1-r_{m,j})}-\sum_{m=1}^b c_{m,j}\right.\nonumber\\ & \quad+\mu_0\alpha_j \left(\sum_{m=1}^b \frac{f_m[G(\omega_m^{{-}1})-1]\omega_m}{\mu_0(1-\omega_m)}\right.\nonumber\\ & \quad\left.\left.+\sum_{l=1}^h\sum_{m=1}^b\frac{c_{m,l}\theta_l[G(r_{m,l}^{{-}1})-1]r_{m,l}}{[(\mu_0-\mu_l) (1-r_{m,l})-\theta_l][\theta_l+\mu_l(1-r_{m,l})]}\right)\right],\quad j=1,\ldots,h. \end{align}

As for the stationary probability $P_0$, it may be determined from the normalization equation

(53)\begin{align} P_0& = 1-\sum_{n=1}^\infty P_n-\sum_{n=1}^\infty\sum_{j=1}^hQ_{n,j}-\sum_{j=1}^hQ_{0,j}\nonumber\\ & =1-\sum_{j=1}^h\sum_{m=1}^b\frac{c_{m,j}[G(r_{m,j}^{{-}1})-1](\mu_0-\mu_j)r_{m,j}}{[(\mu_0-\mu_j) (1-r_{m,j})-\theta_j][\theta_j+\mu_j(1-r_{m,j})]}-\sum_{m=1}^b\frac{f_m[G(\omega_m^{{-}1})-1] \omega_m}{\mu_0(1-\omega_m)^2}\nonumber\\ & \quad -\sum_{j=1}^h\frac{1}{\theta_j}\left[\mu_j\sum_{m=1}^b\frac{c_{m,j}[G(r_{m,j}^{{-}1})-1] r_{m,j}}{\theta_j+\mu_j(1-r_{m,j})}-\sum_{m=1}^b c_{m,j}+\mu_0\alpha_j\left(\sum_{m=1}^b \frac{f_m[G(\omega_m^{{-}1})-1]\omega_m}{\mu_0(1-\omega_m)}\right.\right.\nonumber\\ & \quad \left.\left.+\sum_{j=1}^h\sum_{m=1}^b\frac{c_{m,j}\theta_j[G(r_{m,j}^{{-}1})-1]r_{m,j}}{[(\mu_0-\mu_j) (1-r_{m,j})-\theta_j][\theta_j+\mu_j(1-r_{m,j})]}\right)\right]. \end{align}

It is not difficult to see that the most critical step in the above queueing analysis is to construct a set of linear algebraic equations based on a discrete dynamic system and its roots of the characteristic equation. After solving the system of linear equations, we may combine Eqs. (48) and (49) to compute the sojourn time of an arbitrary customer. Moreover, the expected queue length can also be numerically obtained to verify that whether the expected queue length and the average sojourn time satisfy Little's Law.

5. Sojourn time distribution in the system

Here, the main purpose of studying the sojourn time of an arbitrary customer is to provide a way to verify the correctness of the algorithm for computing the queue-length distribution.

The time that a customer spends in the system, from the instant of its arrival to the queue to the instant of its departure from the server, is called the sojourn time. Denote the random variable that describes the quantity mentioned above by $W_R$, and its mean value by ${\rm E}[W_R]$. At the time of an arbitrary test customer's arrival, it sees all the customers that were already in the system plus all other customers in front of it arriving in the same batch. Let $g^-_k$ ($k=0,1,\ldots,b-1$) be the probability of $k$ number of customers ahead of a randomly selected test customer within the batch. Following Burke [6], we have $g^-_k=({1}/{\bar g})\sum _{\tau =k+1}^{\infty }g_\tau$. Next, we obtain the probability distribution function of $W_R$ by conditioning on the fact that at a batch arrival epoch, the server is serving in vacation mode or in normal mode. Specifically, we need to consider two cases:

(1) Suppose the number of customers that arrive in the same bulk as the test customer, but enter service before the test customer is $k$ $(k=0,1,\ldots,b-1)$. Additionally, we further assume that the test customer finds $n$ $(n\ge 0)$ customers already in the system, and the server is on type $j$ vacation upon its arrival. Thus, if $n+k+1$ customers are served before the single ongoing vacation ends, then the test customer's sojourn time is the sum of $n+k+1$ independent exponential service times with common mean $1/\mu _j$. On the other hand, if $i$ $(i=0,1,\ldots,n+k)$ customers are served before the single ongoing vacation ends, and the rest are served in a normal busy period, then the sojourn time of the test customer is the sum of the remaining vacation time and $n+k+1-i$ exponential service times with intensity $\mu _0$.

(2) Suppose a randomly selected test customer's position in an arrival group is $k+1$ $(k=0,1,\ldots,b-1)$, and it finds the server has got back to normal service mode. Further, if there are $n$ customers already in the system when the test customer arrives, then the test customer's sojourn time is equal to the sum of $n+k+1$ regular service times. Let $S_{v,j}$ denote the service time of the $v$th customer in the type $j$ vacation, and $\overline V_j$ represents the remaining vacation time of type $j$. Taking into account every scenario mentioned above, we can conclude

(54)\begin{align} & W_R(t)=\Pr\{W_R\le t\}\nonumber\\ & \quad = \sum_{j=1}^h\sum_{n=0}^\infty Q_{n,j}^-\sum_{k=0}^{b-1}g_k^-\left[\Pr\left\{\left.\sum_{v=1}^{n+k+1}S_{v,j}\le t\, \right|\overline{V}_j\ge\sum_{v=1}^{n+k+1}S_{v,j}\right\}\Pr\left\{\overline{V}_j\ge\sum_{v=1}^{n+k+1}S_{v,j}\right\}\right. \nonumber\\ & \qquad \left.+\sum_{i=0}^{n+k}\Pr\left\{\left.\overline{V}_j+ \sum_{v=i+1}^{n+k+1}S_{v,0}\le t\,\right|\sum_{v=1}^iS_{v,j}\le\overline{V}_j< \sum_{v=1}^{i+1} S_{v,j}\right\}\Pr\left\{\sum_{v=1}^iS_{v,j}\le\overline{V}_j< \sum_{v=1}^{i+1} S_{v,j}\right\}\right]\nonumber\\ & \qquad +\sum_{n=0}^\infty P_n^-\sum_{k=0}^{b-1}g_k^-\Pr\left\{\sum_{v=1}^{n+k+1}S_{v,0}\le t\right\}\nonumber\\ & \quad =\sum_{j=1}^h\sum_{n=0}^\infty Q_{n,j}^-\sum_{k=0}^{b-1}g_k^-\left[\int_0^t\frac{\mu_j(\mu_jx)^{n+k}}{(n+k)!}\,e^{-(\mu_j+ \theta)x}\,{\rm d}x\vphantom{\left[\sum_{v=0}^{n+k-i}\right]}\right.\nonumber\\ & \qquad \left.+\sum_{i=0}^{n+k}\int_0^t\theta_j\,e^{-(\theta_j+\mu _j)x}\frac{(\mu_jx)^i}{i!}\left[1-e^{-\mu_0(t-x)}\sum_{v=0}^{n+k-i}\frac{(\mu_0(t-x))^v}{v!}\right]{\rm d}x\right]\nonumber\\ & \qquad +\sum_{n=0}^\infty P_n^-\sum_{k=0}^{b-1}g_k^-\int_0^t\frac{\mu_0(\mu_0x)^{n+k}}{(n+k)!}\,e^{-\mu_0x}\,{\rm d}x. \end{align}

Let $W_R^*(s)=\int _0^\infty e^{-st}{\rm d}W_R(t)$ be the L.S.T. of $W_R(t)$. From the convolution property of the Laplace transform, we have

(55)\begin{align} W_R^*(s)& = \sum_{j=1}^h\sum_{n=0}^\infty Q_{n,j}^-\sum_{k=0}^{b-1}g_k^-\left[\left( \frac{\mu_j}{s+\mu_j+\theta_j}\right)^{n+k+1}+\sum_{i=0}^{n+k} \frac{\theta_j\mu_j^i}{(s+\mu_j+\theta_j)^{i+1}} \left(\frac{\mu_0}{s+\mu_0}\right)^{n+k+1-i}\right]\nonumber\\ & \quad +\sum_{n=0}^\infty P_n^-\sum_{k=0}^{b-1}g_k^-\left(\frac{\mu_0}{s+\mu_0}\right)^{n+k+1}. \end{align}

By differentiation Eq. (55) with respect to $s$, and setting $s=0$, the expectation of $W_R$ is given by

(56)\begin{align} {\rm E}[W_R]& = \sum_{j=1}^h\sum_{n=0}^\infty Q_{n,j}^-\sum_{k=0}^{b-1}g_k^-\left[\frac{(n+k+1)\mu_j^{n+k+ 1}}{(\mu_j+\theta_j)^{n+k+2}}+ \sum_{i=0}^{n+k}\left(\frac{(i+1)\mu_j^i\theta_j}{(\mu_j+\theta_j)^{i+2}}+ \frac{\mu_j^i\theta_j(n+k+1-i)}{\mu_0(\mu_j+\theta_j)^{i+1}}\right)\right]\nonumber\\ & \quad +\sum_{n=0}^\infty P_n^-\sum_{k=0}^{b-1}g_k^-\frac{n+k+1}{\mu_0}. \end{align}

It is well known that after obtaining the average queue length $L_s$, the mean sojourn time of an arbitrary customer can be evaluated by Little's Law, ${\rm E}[W_R]=L_s/\lambda \overline {g}$. However, Eq. (56) provides us with another semi-analytical way to get this performance measure. Therefore, we say that comparing the consistency of the calculation results of Eq. (56) and Little's Law is an effective means to test the feasibility of this method.

6. Numerical illustrations

To summarize, the procedure for the calculation of the queue-length distribution adopted here consists of the following three-step algorithm.

Step 1. Set the values for the parameters of the system based on the stability condition, and find the roots of the characteristic equations (23) and (35), respectively.

Step 2. Combining Eqs. (29), (46), and (47) gives a system of $bh+b$ linear equations in $bh+b$ unknowns. Since Eq. (29) does not contain the unknown vector $\bf f$, we have to expand the coefficient matrix of Eq. (29) as follows

$$\left(\begin{array}{ccccccc} \bf R_1 & \bf 0 & \bf 0 & \cdots & \bf 0 & \bf 0 & \bf 0\\ \bf 0 & \bf R_2 & \bf 0 & \cdots & \bf 0 & \bf 0 & \bf 0\\ \bf 0 & \bf 0 & \bf R_3 & \cdots & \bf 0 & \bf 0 & \bf 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ \bf 0 & \bf 0 & \bf 0 & \cdots & \bf R_{h-1} & \bf 0 & \bf 0\\ \bf 0 & \bf 0 & \bf 0 & \cdots & \bf 0 & \bf R_h & \bf 0\\ \end{array}\right).$$

Step 3. Solve the above system of linear equations and obtain the stationary queue-length distributions at different epochs by using Eqs. (48), (49), (50), (51), (52), and (53).

The following numerical examples are presented to illustrate the application of our analysis. Consider a GI$^X$/M/1 queue in which the exponential service distribution in a normal busy period has mean service time equal to $0.2$ (i.e., $\mu _0=5$), and in which the server's vacation time has a hyper-exponential distribution function $V(t)=0.8(1-e^{-0.1t})+0.2(1-e^{-0.25t})$. If such a distribution is used to model a server's vacation, then a server entering a maintenance phase will, with probability $80\%$, receive imperfect preventive repair that is exponentially distributed with parameter $\theta _1=0.1$ and then exit the type 1 vacation mode, or else, with probability $20\%$ receive a preventive replacement of a specific component that is exponentially distributed with parameter $\theta _2=0.25$ and then exit the type 2 vacation mode. Furthermore, the exponential service time distributions have an average service time of $0.5$ in type 1 vacation, and $0.8$ in type 2 vacation, i.e., $\mu _1=2$ and $\mu _2=1.25$. For the batch arrival process, we consider four different cases listed below:

Case 1. Investigate a batch arrival Poisson queue. We fix the value of arrival rate $\lambda =0.6$, and illustrate a numerical example by assuming that the p.m.f. of batch size $X$ has the constant value $1/12$ for all possible values $k$ of $X$, $k=1,2,\ldots,12$. This leads to $\overline {g}=6.5$, $\rho =0.78$ and $a^*(s)={0.6}/{(s+0.6)}$.

Case 2. The numerical results in this case were obtained by assuming that the group size distribution is normalized Poisson with p.m.f. $g_k={e^{-0.8}(0.8)^k}/\left({k!\sum _{n=1}^{13} {e^{-0.8}(0.8)^n}/{n!}}\right)$, $k=1,2,\ldots,13$, and the inter-batch arrival time follows a phase-type distribution with an irreducible representation $(\bf \sigma, \bf L)$ of order two, where $\bf \sigma =(0.7, 0.3)$, $\bf L=\left (\begin {smallmatrix} -3 & 2.7 \\ 4.5 & -21 \end {smallmatrix}\right )$. This leads to $\lambda =2.699045$, $\overline {g}=1.452773$, $\rho =0.784220$, and $a^*(s)={(516s+5085)}/{5(20s^2+480s+1017)}$.

Case 3. Suppose that the inter-batch arrival time distribution obeys a deterministic distribution with mean $1.25$, and the arriving batch size follows a 1–3–6–9 distribution with p.m.f. $g_1=0.1$, $g_3=0.25$, $g_6=0.45$, $g_9=0.2$. This leads to $\lambda =0.8$, $\overline {g}=5.35$, $\rho =0.856$, and $a^*(s)=e^{-1.25s}$.

Case 4. Model the batch arrival process by inverse Gaussian distribution with p.d.f. $a(t)=({0.75}/{\sqrt {2\pi t^3}})e^{-{(t-0.75)^2}/{2t}}$, $t>0$, and for which the p.m.f. of bulk size is $g_k={0.35(1-0.35)^{k-1}}/{(1-(1-0.35)^{10})}$, $k=1,\ldots,10$. This leads to $\lambda =4/3$, $\overline {g}=2.720678$, $\rho =0.725514$, and $a^*(s)=e^{0.75-0.75\sqrt {1+2s}}$.

We have found that in the above cases not all of the distribution functions for the inter-batch arrival time have rational Laplace–Stieltjes transforms. This makes it hard or even impossible to use mathematical software to get the roots of the characteristic equations. To overcome this difficulty, the classical technique of Pad$\acute {\rm e}$ approximation is employed to approximate a given Laplace–Stieltjes transform by a suitable rational function $\mathcal {R}(s)$. In the numerical experiments, the MATHEMATICA build-in command “PadeApproximant” can be used to easily generate the Pad$\acute {\rm e}$ approximants. Thus, the Laplace–Stieltjes transforms of the inter-batch arrival time distributions in Case 3 and Case 4 can be approximated as

(57)\begin{align} & e^{{-}1.25s}\approx\frac{\begin{array}{c} 1-0.588235s+0.160846s^2-0.026808s^3+0.002992s^4-0.000230s^5\\ + 0.000012s^6-3.89184\times10^{{-}7}s^7 +6.08099\times10^{{-}9}s^8\end{array}}{\begin{array}{c} 1+0.661765s+0.206801s^2+0.040211s^3+0.005385s^4+0.000518s^5\\ + 0.000036s^6+1.75133\times10^{{-}6}s^7+5.4729\times10^{{-}8}s^8+8.44583\times10^{{-}10}s^9 \end{array}}, \end{align}

(58)\begin{align} & e^{0.75-0.75\sqrt{1+2s}}\approx\frac{1+4.95748s+8.70859s^2+6.32864s^3+1.57069s^4}{\begin{array}{c} 1+5.70748s+12.333s^2+12.5594 s^3+6.1046s^4+1.27288s^5\\ + 0.092802s^6+0.003849s^7 \end{array}}. \end{align}

As for the issues on choosing the best Padé approximant, readers may refer to the work done by Singh et al. [24]. Substituting the values of parameters into Eqs. (23) and (35), and using the Padé approximants for $a^*(s)$, we may successfully obtain all roots of the characteristic equations inside the unit circle. For a visual illustration, we plot these roots in Figure 1(a–d). We may observe that the imaginary roots of the characteristic equations always exist in complex conjugate pairs. If $b$ is an odd number, there are ${(b-1)}/{2}$ pairs of complex conjugate roots for each characteristic equation. In contrast, if $b$ is an even number, there are ${(b-2)}/{2}$ pairs of complex conjugate roots for each characteristic equation. Using these roots, we establish the coefficient matrix of the system of $bh+b$ linear equations that determines the unknown constants $c_{m,j}$ and $f_m$.

Figure 1. All the roots of the characteristic equations (23) and (35) lie inside the unit circle for different cases. (a) Exponential inter-batch arrival time and uniformly distributed batch size. (b) PH inter-batch arrival time and normalized Poisson distributed batch size. (c) Deterministic inter-batch arrival time and 1–3–6–9 distributed batch size. (d) Inverse Gaussian inter-batch arrival time and normalized geometrically distributed batch size.

Once $c_{m,j}$ and $f_m$ are obtained, the queue-length distributions can be computed by some basic algebraic operations. Due to lack of space, the queue-length distributions at different epochs presented in Tables 1–4 are reported to six decimal places. The main objective of presenting such results in tabular form is to show that the semi-analytical method has good numerical tractability. Furthermore, we know an important property of the batch Poisson arrival process is that the distribution of customers seen by a batch arrival to a queueing facility is, stochastically, the same as the limiting distribution of customers at that facility (i.e., PASTA property, see Gross and Harris [13]). In other words, once the queueing system has reached a steady-state, the distribution of customers at batch arrival instants is also the same as the distribution of customers at any time instant. There is no doubt that the data in Table 1 confirm this fact. We may also use this property as one of the effective ways to check the accuracy of numerical solutions. Additionally, at the bottom of each table, we show the mean sojourn time of an arbitrary customer calculated by Eq. (56) and Little's Law, respectively. We may find that ${\rm E}[W_R]$ evaluated through Eq. (56) exactly matches the one obtained from Little's Law, namely ${\rm E}[W_R]={\rm E}[W_R]_{\rm Little}$.

Table 1. Queue-length distributions at two different epochs for M$^X$/M/1 queue with single hyper-exponential working vacation.

${\rm E}[W_R]=8.897022$, ${\rm E}[W_R]_{\rm Little}=8.897022$.

Table 2. Queue-length distributions at two different epochs for PH$^X$/M/1 queue with single hyper-exponential working vacation.

${\rm E}[W_R]=6.215603$, ${\rm E}[W_R]_{\rm Little}=6.215603$.

Table 3. Queue-length distributions at two different epochs for D$^X$/M/1 queue with single hyper-exponential working vacation.

${\rm E}[W_R]=6.941437$, ${\rm E}[W_R]_{\rm Little}=6.941437$.

Table 4. Queue-length distributions at two different epochs for Inverse Gaussian$^X$/M/1 queue with single hyper-exponential working vacation.

${\rm E}[W_R]=6.538318$, ${\rm E}[W_R]_{\rm Little}=6.538318$.

7. Conclusions

This paper has studied a more realistic working vacation queue that has never been considered in queueing literature. Using the supplementary variable technique and the shift operator method in the theory of difference equations, we obtain the queue-length distributions at the pre-arrival and arbitrary epochs simultaneously in terms of roots of the associated characteristic equation. Since the root-finding procedure is no more difficult with the help of MATHEMATICA software package, the core of computing the queue-length distribution is to solve a nonhomogeneous system of linear equations in our algorithm. We may see that characteristic roots serve as a bridge in our analysis. During the development of queueing theory, several authors have positively evaluated the advantages of the methods based on the characteristic roots. In comparison with the matrix-geometric method, Daigle and Lucantoni [10] stated that whenever the roots method works, it works blindingly fast and is insensitive to traffic intensity $\rho$, whereas the matrix-geometric method is not. Gouweleeuw [12] also stated that using the roots method to obtain the numerical results for the stationary queue-length probabilities will be more efficient. In fact, the method presented in this paper is straightforward in terms of analysis, notation, and computation. Additionally, this method can effectively avoid discussing the transition probability matrix for the embedded Markov chain, so we do not need to focus our attention on the minimal nonnegative solution to a nonlinear matrix equation arising in many queueing problems. Further, through the extensive numerical experiments, we also see that this algorithm can address many different arrival patterns and is no longer limited to the class of phase-type distributions that have rational Laplace–Stieltjes transform. Finally, incorporating a randomize working vacation policy into a renewal input batch arrival queue is worthy of investigation in our future research.