Stationary analysis of a tandem queue with coupled processors subject to global breakdowns

TL;DR

This paper analyzes a tandem queue with coupled processors under global breakdowns, using Markov modulated random walks and generating functions to derive stationary distributions and system performance insights.

Contribution

It introduces a novel analysis of a tandem queue with coupled processors and global breakdowns, employing power series and boundary value problem methods.

Findings

Power series expansions of stationary queue lengths obtained.

Generating function derived via Riemann-Hilbert boundary value problem.

Numerical results provide system performance insights.

Abstract

We consider a tandem queue with coupled processors, which is subject to global breakdowns. When the network is in the operating mode and both queues are non empty, the total service capacity is shared among the stations according to fixed proportions. When one of the stations becomes empty, the total service capacity is given to the non-empty station. Moreover, arrival rates depend on the state of the network. The system is described by a Markov modulated random walk in the quarter plane representing the number of jobs in the two stations and the state of the network. By applying the generating function approach, we first apply the power series approximation method to obtain power series expansions of the generating function of the stationary queue lengths for both network states. Then, we also provide a way to derive the generating function of the stationary queue lengths for both network states in terms of the solution of a Riemann-Hilbert boundary value problem. Numerical results are obtained to show insights in the system performance.

Figures

Equations

$$\begin{array}[]{rl}(\lambda_{0}+\gamma)\pi_{0}(0,0)=&\nu_{2}\pi_{0}(0,1)+\tau\pi_{1}(0,0),\\ (\lambda_{0}+\gamma+\nu_{2})\pi_{0}(0,1)=&\nu_{1}\pi_{0}(1,0)+\nu_{2}\pi_{0}(0,2)+\tau\pi_{1}(0,1),\\ (\lambda_{0}+\gamma+\nu_{2})\pi_{0}(0,k)=&p\nu_{1}\pi_{0}(1,k-1)+\nu_{2}\pi_{0}(0,k+1)+\tau\pi_{1}(0,k),\,k\geq 2\\ \end{array}$$

$$\begin{array}[]{l}(\lambda_{0}+\gamma+\nu_{1})\pi_{0}(n,0)=\lambda_{0}\pi_{0}(n-1,0)+(1-p)\nu_{2}\pi_{0}(n,1)+\tau\pi_{1}(n,0),\vspace{2mm}\\ (\lambda_{0}+\gamma+p\nu_{1}+(1-p)\nu_{2})\pi_{0}(n,1)=\lambda_{0}\pi_{0}(n-1,1)+\nu_{1}\pi_{0}(n+1,0)\\ +(1-p)\nu_{2}\pi_{0}(n,2)+\tau\pi_{1}(n,1),\,n\geq 1,\vspace{2mm}\\ (\lambda_{0}+\gamma+p\nu_{1}+(1-p)\nu_{2})\pi_{0}(n,k)=\lambda_{0}\pi_{0}(n-1,k)+p\nu_{1}\pi_{0}(n+1,k-1)\\ +(1-p)\nu_{2}\pi_{0}(n,k+1)+\tau\pi_{1}(n,k),\,n\geq 1,k\geq 2\end{array}$$

$$\begin{array}[]{c}(\lambda_{1}+\tau)\pi_{1}(n,k)=\gamma\pi_{0}(n,k)+\lambda_{1}\pi_{1}(n-1,k).\end{array}$$

$${\Pi }_i(x,y)=\sum _{n=0}^{\infty }\sum _{k=0}^{\infty }{\pi }_i(n,k)x^n{y}^k,i=0,1,|x|\le 1,|y|\le 1.$$

$$\begin{array}[]{rl}R(x,y)\Pi_{0}(x,y)=&A(x,y)\Pi_{0}(x,0)+B(x,y)\Pi_{0}(0,y)\\ &+C(x,y)\Pi_{0}(0,0)+\tau xy\Pi_{1}(x,y),\vspace{2mm}\\ \ \Pi_{1}(x,y)=&\frac{\gamma}{D(x)}\Pi_{0}(x,y),\end{array}$$

$$\begin{array}[]{rl}R(x,y)=&xy(\lambda_{0}(1-x)+\gamma)+\nu_{1}py(x-y)+\nu_{2}(1-p)x(y-1),\\ A(x,y)=&(1-p)[\nu_{2}x(y-1)+\nu_{1}y(y-x)],\\ B(x,y)=&-\frac{p}{1-p}A(x,y),\\ C(x,y)=&\nu_{1}(1-p)y(x-y)+\nu_{2}px(y-1).\end{array}$$

$$\begin{array}[]{c}\Pi_{0}(x,y)[D(x)R(x,y)-\tau\gamma xy]=D(x)\{A(x,y)\Pi_{0}(x,0)+B(x,y)\Pi_{0}(0,y)\\ +C(x,y)\Pi_{0}(0,0)\}.\end{array}$$

$${\Pi }_0(1,1)=\frac{\tau }{\tau +\gamma },{\Pi }_1(1,1)=\frac{\gamma }{\tau +\gamma }.$$

$${\Pi }_0(s(y),y)=\frac{D(s(y))C(s(y),y)}{D(s(y))R(s(y),y)-\tau \gamma s(y)y}{\Pi }_0(0,0).$$

$${\Pi }_0(0,0)=\frac{\tau }{\tau +\gamma }-\left({\lambda }_0\frac{\tau }{\tau +\gamma }+{\lambda }_1\frac{\gamma }{\tau +\gamma }\right)\left(\frac{1}{{\nu }_1}+\frac{1}{{\nu }_2}\right).$$

$${\lambda }_0\left(\frac{1}{{\nu }_1}+\frac{1}{{\nu }_2}\right)\frac{\tau }{\tau +\gamma }+{\lambda }_1\left(\frac{1}{{\nu }_1}+\frac{1}{{\nu }_2}\right)\frac{\gamma }{\tau +\gamma }<\frac{\tau }{\tau +\gamma },$$

$${\Pi }_0(0,0)=\frac{\tau }{\tau +\gamma }(1-\frac{{\rho }_0\tau +{\rho }_1\gamma }{\tau }).$$

$$y:=\xi (x)=\frac{{\nu }_{2}xD(x)}{xD(x)({\nu }_2+{\lambda }_0(1-x))+{\lambda }_1\gamma x(1-x)},$$

$${\Pi }_0(x,0)=\left(\frac{\tau }{\tau +\gamma }\right)\frac{(1-\frac{{\rho }_0\tau +{\rho }_1\gamma }{\tau }){\nu }_1\xi (x)(\xi (x)-x)}{{\nu }_{2}x(\xi (x)-1)+{\nu }_1\xi (x)(x-\xi (x))}.$$

$$\begin{array}[]{rl}\Pi_{0}(x,y)=&\frac{(\frac{\tau}{\tau+\gamma})(1-\frac{\rho_{0}\tau+\rho_{1}\gamma}{\tau})D(x)\nu_{1}\nu_{2}}{D(x)[\lambda_{0}y(1-x)+\nu_{2}(y-1)]+\lambda_{1}\gamma y(1-x)}\vspace{2mm}\\ &\times\left\{\frac{\xi(x)(y-1)(\xi(x)-x)+y(x-y)\nu_{2}(\xi(x)-1)}{\nu_{2}x(\xi(x)-1)+\nu_{1}\xi(x)(x-\xi(x))}\right\},\vspace{2mm}\\ \Pi_{1}(x,y)=&\frac{\gamma}{D(x)}\Pi_{0}(x,y).\end{array}$$

$$\begin{array}[]{l}[D(x)(\lambda_{0}x(1-x)+\nu_{1}(x-y))+\lambda_{1}\gamma x(1-x)]y\Pi_{0}(x,y)=\vspace{2mm}\\ D(x)\{\Pi_{0}(0,y)[\nu_{2}x(1-y)+\nu_{1}y(x-y)]+\Pi_{0}(0,0)\nu_{2}x(y-1)\}.\end{array}$$

$${\Pi }_0(0,0)={\Pi }_0(0,y)[1-\frac{{\nu }_{1}y(u(y)-y)}{{\nu }_{2}u(y)(1-y)}].$$

$$\begin{array}[]{rl}\Pi_{0}(x,y)=&\frac{D(x)\nu_{1}y(x-u(y))}{u(y)[D(x)(\lambda_{0}x(1-x)+\nu_{1}(x-y))+\lambda_{1}\gamma x(1-x)]}\Pi_{0}(0,y),\vspace{2mm}\\ \Pi_{1}(x,y)=&\frac{\gamma}{D(x)}\Pi_{0}(x,y).\end{array}$$

$${\Pi }_j(x,y)=\sum _{m=0}^{\infty }{V}_m^{(j)}(x,y)p^m,j=0,1.$$

$$\begin{array}[]{c}G(x,y)\Pi_{0}(x,y)-G_{10}(x,y)\Pi_{0}(x,0)-G_{00}(x,y)\Pi_{0}(0,0)\vspace{2mm}\\ =pG_{10}(x,y)[\Pi_{0}(x,y)-\Pi_{0}(x,0)-\Pi_{0}(0,y)+\Pi_{0}(0,0)],\end{array}$$

$$\begin{array}[]{rl}G(x,y)=&D(x)[\lambda_{0}y(1-x)+\nu_{2}(y-1)]+\lambda_{1}\gamma y(1-x),\\ G_{10}(x,y)=&D(x)[\nu_{2}(y-1)-\nu_{1}y(1-yx^{-1})],\\ G_{00}(x,y)=&D(x)\nu_{1}y(1-yx^{-1}).\end{array}$$

$$\begin{array}[]{rl}V_{0}^{(0)}(x,y)=&\frac{D(x)\nu_{1}\nu_{2}(\tilde{Y}(x)-x)(\tilde{Y}(x)(y-1)+x-y)}{G(x,y)[\nu_{2}x(\tilde{Y}(x)-1)-\nu_{1}\tilde{Y}(x)(x-\tilde{Y}(x))]}V^{(0)}(0,0),\vspace{2mm}\\ V_{m}^{(0)}(x,y)=&\frac{G_{10}(x,y)}{G(x,y)}Q_{m-1}(x,y),\,m>0,\vspace{2mm}\\ V_{m}^{(1)}(x,y)=&\frac{\gamma}{D(x)}V_{m}^{(0)}(x,y),\,m\geq 0,\end{array}$$

$$\begin{array}[]{rl}\tilde{Y}(x)=&\frac{\nu_{2}D(x)}{D(x)(\nu_{2}+\lambda_{0}(1-x))+\lambda_{1}\gamma(1-x)},\\ Q_{m}(x,y)=&V_{m}^{(0)}(x,y)-V_{m}^{(0)}(x,\tilde{Y}(x))-V_{m}^{(0)}(0,y)+V_{m}^{(0)}(0,\tilde{Y}(x)),\,m\geq 0,\end{array}$$

$$\begin{array}[]{rr}V_{m}^{(0)}(x,y)G(x,y)=&G_{10}(x,y)[V_{m}^{(0)}(x,0)+P_{m-1}(x,y)]\\ &+G_{00}(x,y)V_{m}^{(0)}(0,0),\,m\geq 0,\end{array}$$

$$P_m(x,y)=V_m^{(0)}(x,y)-V_m^{(0)}(x,0)-V_m^{(0)}(0,y)+V_m^{(0)}(0,0),$$

$$\tilde{Y}(x)=\frac{{\nu }_{2}D(x)}{D(x)({\nu }_2+{\lambda }_0(1-x))+{\lambda }_1\gamma (1-x)}.$$

$$V_m^{(0)}(x,0)=-\frac{G_{00}(x,\tilde{Y}(x))}{G_{10}(x,\tilde{Y}(x))}{V}_m^{(0)}(0,0)-P_{m-1}(x,y).$$

$$\begin{array}[]{rl}V_{0}^{(0)}(0,0)=&\frac{\tau}{\tau+\gamma}-\left(\lambda_{0}\frac{\tau}{\tau+\gamma}+\lambda_{1}\frac{\gamma}{\tau+\gamma}\right)\left(\frac{1}{\nu_{1}}+\frac{1}{\nu_{2}}\right),\\ V_{m}^{(0)}(0,0)=&0,\,m>0.\end{array}$$

$$\begin{array}[]{rl}E(Q_{1})=&\sum_{m=0}^{\infty}p^{m}\frac{\partial}{\partial x}[V_{m}^{(0)}(x,1)+V_{m}^{(1)}(x,1)]|_{x=1}\\ =&\frac{\lambda_{1}\gamma}{\tau(\tau+\gamma)}+(1+\frac{\gamma}{\tau})\sum_{m=0}^{\infty}p^{m}\frac{\partial}{\partial x}V_{m}^{(0)}(x,1)|_{x=1},\vspace{2mm}\\ E(Q_{2})=&\sum_{m=0}^{\infty}p^{m}\frac{\partial}{\partial y}[V_{m}^{(0)}(1,y)+V_{m}^{(1)}(1,y)]|_{y=1}\vspace{2mm}\\ =&(1+\frac{\gamma}{\tau})\sum_{m=0}^{\infty}p^{m}\frac{\partial}{\partial y}V_{m}^{(0)}(1,y)|_{y=1},\end{array}$$

$$v_{m,1}=\frac{\partial }{\partial x}{V}_m^{(0)}(x,1){|}_{x=1},v_{m,2}=\frac{\partial }{\partial y}{V}_m^{(0)}(1,y){|}_{y=1}.$$

Full text

1 Introduction

Queueing networks with service interruptions are known to be adequate models to handle realistic problems in manufacturing, telecommunications etc. Despite their great importance, there has been done very few works involving more than one queue, since even under favorable assumptions, the existence of service interruptions destroys the separability [5] of the appropriate multidimensional Markov process and renders its solution intractable.

The problem becomes even more challenging when we further assume that the nodes are interacting with each other. With the term interaction", we mean that the service rate at a node depends on the state of the other nodes, i.e., a networks with coupled processors [[20](#bib.bib20), [31](#bib.bib31), [26](#bib.bib26)]. For such systems, which do not possess the product form" solution, analytic methods have been developed in [21, 10].

In this work, we go one step further and consider a two-node tandem network with coupled processors, which is subject to global (i.e., network) breakdowns. It is assumed that jobs arrive at the first station according to a Poisson process according to the state of the network, and require service at both stations before leaving the system. The amounts of work that a job requires at each of the stations are independent, exponentially distributed random variables. When the network is in the operating mode, and both stations are nonempty, the total service capacity is shared between the stations according to fixed proportions.

When one of the stations becomes empty, the total service capacity is given to the non-empty station. When the network is in the setup mode after a failure occurrence, both stations stop working111This is natural when we are dealing with queues in series., for an exponentially distributed time period. During setup period, jobs continue to arrive, but now at a decreased rate in order to avoid further congestion222Such an operation can be performed by a central scheduler, responsible for the congestion management.. Such a system is fully described by a random walk in the quarter plane (representing the number of jobs in each node), which is modulated by a two state Markov chain (representing the state of the network).

For such a network we provide to different approaches, named the Power Series Approximation (PSA) method, and the theory of Riemann-Hilbert boundary value problems (BVP) to investigate its stationary behavior.

1.1 Related work

Most of the existing studies involving breakdowns have concentrated on models with a single job queue served by one or more processors, e.g. [32, 28]. Other related results where jobs from a single source are directed to one of several parallel queues, and where breakdowns result in the loss of jobs, or the direct transfer to other queues are given in [33, 27, 34]. A network of two nodes subject to breakdowns and repairs was also analysed in [29]. Approximate solutions to obtain performance measures, based on replacing an interruptable server with an uninterruptable but slower one, choosing the new service rate without affecting the overall service capacity was given in [30, 37].

Queues with coupled processors were initially studied in [20]. To gain quantitative insights about the queueing process, they derived a solution for the generating function of the stationary joint queue-length distribution using the theory of Riemann-Hilbert boundary value problems. Later, in [10] a systematic and detailed study of the technique of reducing a two dimensional functional equation of a random walk or queueing model to a boundary value problem was presented, while several numerical issues were discussed. Important generalizations were given in [9, 4, 24, 26, 31, 14, 13, 16, 17, 18, 19] (not exhaustive list) where various two-dimensional queueing models with the aid of the theory of Riemann (-Hilbert) boundary value problems.

A tandem queue with two coupled processors was analyzed in [31], while later computational issues, as well as asymptotic results were discussed and presented in [26, 23], respectively. There, it was shown that the problem of finding the bivariate generating function of the stationary joint queue length distribution is reduced to a Riemann-Hilbert boundary value problem.

Applications of coupled processor models arise systems where limited resources are dynamically shared among processors, e.g., in data transfer in bidirectional cable and data networks [26], in bandwidth sharing of data flows [23], in the performance modeling of device-to-device communication [38], to model the complex interdependence among transmitters due to interference [7, 18, 16, 19], as well as in assembly lines in manufacturing [2].

Other approaches to analyze two-dimensional queueing models have been developed in [1] (compensation method), and [6, 25] (power series algorithm). In the latter one, power series expansions of steady-state probabilities as functions of a certain parameter of the system (usually the load) were derived. Recently, the authors in [39] studied generalized processor sharing queues and introduced an alternative method, by providing power series expansions of the generating function; see also [36, 35, 15].

1.2 Our contribution

In this work we focus on the stationary analysis of two-node tandem queue with coupled processors and network breakdowns.

Applications

Potential application of our system are found in systems with limited capacity, which must be shared in multiple operations. For example, in the modeling of virus attacks or other malfunctions in cable access networks regulated by a request-grant mechanism.

Another application of the model can be found in manufacturing [2], and in particular in an assembly line. There, two operations on each job must be performed using a limited service capacity. To increase the network throughput, we couple the service rates at each of the operations, and thus we use the service capacity of an operation for which no jobs are waiting for the other operation. Such a system is heavily affected by the presence of failures during the job processing, which in turn will definitely deteriorate the system throughput.

Fundamental contribution

Based on the generating function approach, we apply two different methods to investigate the stationary behaviour of the underline Markov modulated random walk in the quarter plane. First, we apply the power series approximation method, initially introduced in [39] (see also [11, 35, 36]), for two-parallel generalized processor sharing queues, described by a typical random walk in the quarter plane (RWQP). In this work we show that this method is still valid for related Markov modulated RWQP, and thus, extend the class of models that can be applied. Under such a method we obtain power series expansions of the probability generating function (pgf) of the joint stationary distribution for either state of the network. A recursive technique to derive their coefficients is also presented.

Secondly, we also derive the pgfs of the joint stationary queue-length distribution for either state of the network with the aid of the theory of Riemann-Hilbert boundary value problems [20, 21, 10]. It is seen that applying the theory of boundary value problems some further technical difficulties are also arise. More precisely, by applying the generating function approach we first come up with a system of functional equations, which is then reduced to a single fundamental equation.

The rest of the paper is summarized as follows. In Section 2 we present the mathematical model in detail and obtain the functional equations along with some preliminary results. Section 3 is devoted to the analysis of the two extreme cases where the total capacity is allocated to one of the two stations, even if both stations are nonempty. In Section 4 we apply the power series approximation method for the case where the service capacity is shared by the two stations, where in Section 5, we provide a complete analysis on how to obtain the pgfs of the stationary joint queue length distribution in terms of a solution of a Riemann-Hilbert boundary value problem. Numerical validations of the performance metrics obtained by using the PSA and the BVP approaches are given in 6.

2 The model and the functional equations

Consider a two-stage tandem queue, where jobs arrive at queue 1 according to a Poisson process with rate depending on the state of the network. In particular, the network is subject to breakdowns which occur according to a Poisson process with rate $\gamma$. When a breakdown occurs, both stations stop working for an exponentially distributed time period with rate $\tau$. Thus, the network alternates between the operating mode and the setup mode. Denote by $C(t)$ the state of the network at time $t$, with $C(t)=0$ (resp. 1), when network is in operating (resp. setup). When $C(t)=i$, jobs arrive in station 1 according to a Poisson process with rate $\lambda_{i}$, $i=0,1$.

Each job demands service at both queues before departing from the network. More precisely, at station $j$, a job requires an exponentially distributed amount of service with parameter $\nu_{j}$, $j=1,2$. The total service capacity of the tandem network equals one unit of work per time unit. In particular, when both stations are non-empty, station $j$ is served at a rate $\phi_{j}$, $j=1,2$, and at a rate 1 when it is the only non empty. Without loss of generality we assume hereon that $\phi_{1}=p$ and $\phi_{2}=1-p$, where $0\leq p\leq 1$.

Let $Q_{j}(t)$, $j=1,2,$ be the number of customers at queue $j$ at time $t$. Under usual assumptions the stochastic process $X(t)=\{(C(t),Q_{1}(t),Q_{2}(t));t\geq 0\}$ is an irreducible and aperiodic continuous time Markov chain with state space $E=\{0,1\}\times\mathbb{Z}^{+}\times\mathbb{Z}^{+}$. Denote by $\pi_{i}(n,k)$ the stationary probability of having $n$ and $k$ customers at stations 1 and 2, respectively, when the network is in state $i$. The balance equations are given by

[TABLE]

[TABLE]

[TABLE]

Define the probability generating functions of the joint stationary queue length distribution

[TABLE]

Using the balance equations we obtain after some algebra the following system of functional equations

[TABLE]

where, $D(x)=\lambda_{1}(1-x)+\tau$ and

[TABLE]

Our aim is to solve the system of functional equations (4). Substituting the second in (4) to the first one, we obtain the following fundamental functional equation

[TABLE]

Clearly, $\Pi_{0}(1,1)+\Pi_{1}(1,1)=1$, while by using the second in (4) we obtain $\Pi_{1}(1,1)=\frac{\tau}{\gamma}\Pi_{0}(1,1)$. Thus, the probabilities of the network state are easily given by

[TABLE]

Let $x=s(y):=\frac{\nu_{1}y^{2}}{\nu_{1}+\nu_{2}(1-y)}$. Note that $A(s(y),y)=B(s(y),y)=0$. Then, using (6) we obtain

[TABLE]

Letting $y\to 1$ in (7) we obtain

[TABLE]

Note that (8) implies that our network is stable when

[TABLE]

which can be explained by realizing that the left hand side of (9) equals the amount of work brought into the system per time unit, and in order the system to be stable, should be less than the amount of work departing the system per time unit.

Let $\rho_{kj}:=\frac{\lambda_{k}}{\nu_{j}}$, $k=0,1$, $j=1,2$, and $\rho_{k}=\rho_{k1}+\rho_{k2}$, $k=1,2$. Then, (8) is rewritten as

[TABLE]

Remark 1

Note that $\lambda_{0}\left(\frac{1}{\nu_{1}}+\frac{1}{\nu_{2}}\right)\frac{\tau}{\tau+\gamma}$ (resp. $\lambda_{1}\left(\frac{1}{\nu_{1}}+\frac{1}{\nu_{2}}\right)\frac{\gamma}{\tau+\gamma}$) refers to the amount of work that arrive at the system per time unit when the network is in the operating mode (resp. in the setup mode), while a job can depart from the network only when it is in the operating mode, and this is happening with probability $\tau/(\tau+\gamma)$.

3 The cases $p=0$ and 1

When $p=0$ (resp. $p=1$), the model can be seen as a tandem queues served by a single server, in which preemptive priority is given to station 2 (resp. station 1). It is easily seen that in such cases, the functional equation (4) can be easily solved since either the coefficient of $\Pi_{0}(0,y)$ (when $p=0$), or the one of $\Pi_{0}(x,0)$ (when $p=1$) is equal to zero.

In case $p=0$ (i.e., $B(x,y)=0$), upon a service completion in station 1, the server continues serving the customer in station 2, since station 2 has priority. Thus, in such a case our system reduces to an unreliable queueing system, in which the service time consists of two exponential phases with parameters $\nu_{1}$ and $\nu_{2}$, respectively. Note that for $p=0$,

[TABLE]

vanishes the left-hand side of (4), and yields

[TABLE]

Substituting back in (4) yields

[TABLE]

The case $p=1$ is even more interesting and corresponds to an unreliable tandem queue attended by a single server and preemptive priority for the first station. That is, if upon a customer arrival the server is at the second station, it switches immediately to the first one. Moreover, upon a setup completion, after a global breakdown, the server will start serving at the first station if there are customers waiting. And this is the case even if a breakdown occurs when was serving a customer at the second station. To the author's best knowledge, that case has never considered before. For $p=1$, (4) reduces to

[TABLE]

Let $x:=u(y)$ the unique root of $D(x)(\lambda_{0}x(1-x)+\nu_{1}(x-y))+\lambda_{1}\gamma x(1-x)=0$ inside the unit circle. Then, the right-hand side should also vanish and thus,

[TABLE]

Substituting back in (11) yields

[TABLE]

4 The case $0<p<1$: Power Series Approximation in $p$

In the following, we are going to construct a power series expansion of the pgf $\Pi_{0}(x,y)$ in $p$ starting by (6). Then, having that result we are able to construct power series expansions of $\Pi_{1}(x,y)$ in $p$ using the second in (4). With that in mind, let

[TABLE]

Our aim in the following, is to obtain $V_{m}^{(j)}(x,y)$, $m\geq 0$, $j=0,1$, by employing an approach similar to the one developed in [11, 39, 15]333See Appendix 0.C for the analyticity of $\Pi_{j}(x,y)$ close to $p=0$.. Equation (6) is rewritten as

[TABLE]

where,

[TABLE]

The major difficulty in solving (6) corresponds to the presence of the two unknown boundary functions $\Pi_{0}(x,0)$, $\Pi_{0}(0,y)$. Having in mind that in the left-hand side of (14) there is only one boundary function, we are able to follow the approach in [39, 11, 15]. The next Theorem summarizes the basic result of this section.

Theorem 4.1

Under stability condition (9),

[TABLE]

where

[TABLE]

and $Q_{-1}(x,y):=0$.

Proof

The proof follows the lines in [39, 15]. Note that $\Pi_{0}(x,y)$ is analytic function of $p$ in a neighbourhood of [math]. We start by expressing $\Pi_{0}(x,y)$ in power series expansion of $p$ by using (14), and equating the corresponding powers of $p$ at both sides. This yields

[TABLE]

where

[TABLE]

with $P_{-1}(x,y)=0$. Note that $G(x,y)=0$ has a unique zero $y=\tilde{Y}(x)$ such that

[TABLE]

It is easy to realize that $|\tilde{Y}(x)|<1$, for $|x|=1$. Substituting in (16) we eliminate its left-hand side yielding

[TABLE]

Substituting (17) back in (16) yields the first in (15). Then, using the second equation in (4) we derive the coefficients $V_{m}^{(1)}(x,y)$ in terms of $V_{m}^{(0)}(x,y)$ as given in the second in (15). From (8) it is readily seen that

[TABLE]

4.1 Performance metrics

We focus on the mean queue lengths given by

[TABLE]

Let

[TABLE]

Truncation of the power series in (18) yields,

[TABLE]

Truncation yields accurate approximations for $p$ close to 0. However, we have to note the actual calculation of the expressions in (18) requires the computation of the first derivatives of $V_{m}^{(0)}(x,y)$ for $m\geq 0$. Although we provided an algorithm (see Theorem 4.1) to calculate these coefficients, the calculation of their first derivatives is far from straightforward due to the extensive use of L’Hopital’s rule.

5 The case $0<p<1$: A Riemann-Hilbert Boundary Value Problem

In the following, we proceed with the determination of $\Pi_{j}(x,y)$, $j=0,1,$ $|x|\leq 1$, $|y|\leq 1$ with the aid of the theory of boundary value problems (BVP). In particular, we first obtain $\Pi_{0}(x,y)$ in terms of the solution of a Riemann-Hilbert boundary value problem by using (6), and then, we use the second in (4), to finally derive $\Pi_{1}(x,y)$.

A key step to analyze the functional equation (6) is the careful examination of the algebraic curve defined by the kernel equation,

[TABLE]

It is easily seen that $H(x,y)$ is a polynomial of third degree in $x$, and of second degree in $y$. The study of $H(x,y)=0$ (see Appendix 0.A) allows to continue the unknown functions $\Pi_{0}(x,0)$ , $\Pi_{0}(0,y)$ analytically outside the unit disk, and to reduce their determination to a Dirichlet boundary value problem.

Let $\mathcal{C}_{x}=\{x\in\mathbb{C}:|x|=1\}$, $\mathcal{C}_{y}=\{y\in\mathbb{C}:|y|=1\}$, $\mathcal{D}_{x}=\{x\in\mathbb{C}:|x|\leq 1\}$, $\mathcal{D}_{y}=\{y\in\mathbb{C}:|y|\leq 1\}$, and denote by $\mathcal{U}^{+}$ (resp. $\mathcal{U}^{-}$) the interior (resp. the exterior) domain bounded by the contour $\mathcal{U}$. Then the following lemma provides information about the location of the zeros of the kernel $H(x,y)$.

Lemma 1

If $y\in\mathcal{C}_{y}$ (resp. $x\in\mathcal{C}_{x}$), $H(x,y)=0$ has a unique root, say $X_{0}(y)\in\mathcal{D}_{x}$ (resp. $Y_{0}(x)\in\mathcal{D}_{y}$).

Proof

See Appendix 0.A

By the implicit function theorem, we see that the algebraic function $Y(x)$ (respectively $X(y)$) defined by $H(x,Y(x))=0$ (resp. $H(X(y),y)=0$) is analytic except at branch points. Denote $X_{1}(y)$, $X_{2}(y)$ the other two in $x$, with $|X_{1}(1)|<|X_{2}(1)|$, by $Y_{1}(x)$ the other one in $y$.

Lemma 2

The algebraic function $Y(x)$, defined by $H(x,Y(x))=0$, has six real positive branch points, and two of them, say $x_{1}$, $x_{2}$, are such that $0=x_{1}<x_{2}<1$.

Proof

See Appendix 0.A.

Define the image contour, $\mathcal{L}=Y_{0}([\overrightarrow{\underleftarrow{0,x_{2}}}])$, where $[\overrightarrow{\underleftarrow{u,v}}]$ stands for the contour traversed from $u$ to $v$ along the upper edge of the slit $[u,v]$ and then back to $u$ along the lower edge of the slit. The following lemma shows that the mapping $Y(x)$, $x\in[0,x_{2}]$, gives rise to the smooth and closed contour $\mathcal{L}$.

Remark 2

The study of $H(x,y)=0$ with respect to $x$ is slightly more difficult, but allows to also show that the algebraic function $X(y)$ has also two real and non-negative branch points inside $\mathcal{D}_{y}$, say $0\leq y_{1}<y_{2}<1$. Similarly, for $y\in[y_{1},y_{2}]$, $X(y)$ lies on a closed contour $\mathcal{M}$. Further details are omitted due to space limitations, and mainly due to the fact that our main contribution relies on the use of PSA method.

Lemma 3

For $x\in[0,x_{2}]$, the algebraic function $Y(x)$ lies on a closed contour $\mathcal{L}$, which is symmetric with respect to the real line and such that

[TABLE]

Proof

Follows directly from the fact that $\Delta(x)$ is negative for $x\in(0,x_{2})$ (see Appendix 0.A).

5.1 A boundary value problem for $\Pi_{0}(0,y)$

For $x\in\mathcal{C}_{x}$, $y=Y_{0}(x)$ we obtain

[TABLE]

where $F(x,y)=\frac{A(x,y)}{1-p}=-\frac{B(x,y)}{p}$. For $y\in\mathcal{C}_{y}$, $x=X_{0}(y)$

[TABLE]

Equation (22) implies that

[TABLE]

which is a meromorphic function444Its possible poles are the zeros of $F(X_{0}(y),y)$ in $\mathcal{L}^{+}\cap\mathcal{C}_{y}^{-}$.. For $y\in\mathcal{L}^{+}\cap\mathcal{C}_{y}^{-}$, $|X_{0}(y)|<1$, and thus, we can construct analytic continuation of $\Pi_{0}(0,y)$ for all $y\in\zeta_{y}:=\mathcal{L}^{+}\cap\mathcal{C}_{y}^{-}$555Note that $\zeta_{y}$ is an empty set when $(1-p)\nu_{2}\leq p\nu_{1}$, since in such a case $\mathcal{L}$ lies entirely inside the unit circle..

Moreover, as $Y_{0}(x)$ is analytic in $\mathcal{D}_{x}-[0,x_{2}]$ [21], $\Pi_{0}(0,Y_{0}(x))$ is meromorphic in $\mathcal{D}_{x}-[0,x_{2}]$, and from (21)

[TABLE]

We therefore have the relation (21), not only for $x\in\mathcal{C}_{x}$ but also for $x\in\mathcal{D}_{x}-[0,x_{2}]$, and by continuity, for $x\in\mathcal{D}_{x}-[0,x_{2}]$ too. Since $\Pi_{0}(x,0)$ is real for $x\in[0,x_{2}]$, we obtain

[TABLE]

or equivalently,

[TABLE]

Thus our problem is reduced to the determination of a function, which is regular for $y\in\mathcal{L}^{+}$, continuous in $\mathcal{L}^{+}\cup\mathcal{L}$ satisfying the boundary condition (24). A standard way to solve this Riemann-Hilbert boundary value problem is to conformally transformed it on the unit circle [21, 26] by introducing the conformal mappings $z=\gamma(y):\mathcal{L}^{+}\to\mathcal{C}_{y}^{+}$, and its inverse $y=\gamma_{0}(z):\mathcal{C}_{y}^{+}\to\mathcal{L}^{+}$666See Appendix 0.B for details on the numerical derivation of the conformal mappings.

Then, we have the following problem: Find a function $T(z)=\Pi_{0}(0,\gamma_{0}(z))$ regular for $z\in\mathcal{C}_{z}^{+}$, and continuous for $z\in\mathcal{C}_{z}\cup\mathcal{C}_{z}^{+}$ such that, $Re(iT(z))=c(\gamma_{0}(z))$, $z\in\mathcal{C}$. Its solution (see [22]) is given by

[TABLE]

where $K$ is a constant. Then, (21) can be used to obtain $\Pi_{0}(x,0)$, $x\in\mathcal{C}_{x}$ 777Note that a similar analysis can be also performed to obtain $\Pi_{0}(x,0)$ in terms of another Riemann-Hilbert boundary value problem. For $x\in\mathcal{C}_{x}^{+}$, $\Pi_{0}(x,0)$ can be derived by using the Cauchy's formula, yielding

[TABLE]

where

[TABLE]

Using (6) we obtain $\Pi_{0}(x,y)$. Then, using the second in (4) we obtain $\Pi_{1}(x,y)$ and all unknowns are fully specified.

6 Numerical Results

We now compare the PSA to the exact results derived by the BVP approach and investigate the influence of some parameters on the mean queue lengths.

Figure 1 depicts the approximations (19) as a function of $p$ for increasing values of $M$. Set $\lambda_{0}=1$, $\lambda_{1}=0.5$, $\tau=4$, $\gamma=2$, $\nu_{1}=4$, $\nu_{2}=5$. The horizontal lines ($M=0$) equal the values for the tandem system with priority for the second queue. Figure 1 confirms that the PSA approximations are accurate for $p$ close to 0, and clearly, more terms provide larger regions for $p$ where the accuracy is good.

In Figure 2 ($\lambda_{0}=1$, $\lambda_{1}=0.5$, $M=3$, $\nu_{1}=4$, $\nu_{2}=5$), we can observe that the increase in $\gamma$ will definitely increase $E(Q_{2})$, since the system switches to the setup mode more frequently.

7 Conclusion

In this work, we investigated the stationary behaviour of an unreliable two-node tandem queue with coupled processors, which is described by a markov modulated RWQP. Based on the generating function approach, we applied the PSA method and obtained power series expansions of the pgfs of the stationary joint queue length distributions for each state of the network. With this result we shown the flexibility of the PSA approach to be applied in more complicated models. Moreover, we also obtained the corresponding pgfs with the aid of the theory of Riemann-Hilbert boundary value problems. By truncating the power series, we find good approximations for the expected number of customers especially when $p$ is close to 0, by comparing them with the ``exact" derivations through BVP.

In the future, we plan to expand our results to networks with more than two nodes, where the theory of BVPs cannot be applied, as well as to consider general routing among the nodes. Moreover, it could be also interesting to compare PSA method with other approximation techniques developed so far [34, 8, 30, 37]. Moreover, it would be interesting to derive asymptotic estimates for the occurrence of large queue lengths due to the presence of failures.

Appendix 0.A Analysis of the kernel

Proof of Lemma 1 Let $u(x,y)=\lambda_{0}(1-x)+\nu_{1}p(1-\frac{y}{x})+\nu_{2}(1-p)(1-\frac{1}{y})+\gamma$. Note that

[TABLE]

Using the principal value argument it is seen that the number of zeros of $H(x,y)$ in $\mathcal{D}_{x}$ equals the number of zeros of $D(x)u(x,y)$ in $\mathcal{D}_{x}$, which is equal to one. $\square$

Proof of Lemma 2 The branch points of the two-valued function $Y(x)$ are the zeros of the discriminant $\Delta(x)=x(f(x)+g(x))$ of $H(x,y)=0$, where

[TABLE]

Note $g(x)=0$ if and only if $x=0$, $x=1$, and

[TABLE]

Let $x_{1}^{*}$, $x_{2}^{*}$ the zeros of (26). Then,

[TABLE]

which means that $x_{1}^{*},x_{2}^{*}>1$. Thus $g(x)=0$ has exactly two roots in $[0,1]$. By using Rouche's theorem, we can show that $D(x)=0$, i.e., $x(f(x)+g(x))=0$ has exactly two zeros in $[0,1]$, and one of them equals $x_{1}=0$.

Appendix 0.B On the derivation of conformal mappings

A detailed approach on how we can numerically obtain the conformal mappings is given in [10, 26]. We summarized the basic steps. First, we need to represent $\mathcal{L}$ in polar coordinates, i.e., $\mathcal{L}=\{y:y=\rho(\phi)\exp(i\phi),\phi\in[0,2\pi]\}.$ Since $0\in\mathcal{L}^{+}$, for each $y\in\mathcal{L}$, we can have a relation between its absolute value and its real part, i.e., $|y|^{2}=m(Re(y))$. Given the angle $\phi$ of some point on $\mathcal{L}$, the real part of this point, say $\delta(\phi)$, is the solution of $\delta-\cos(\phi)\sqrt{m(\delta)}$, $\phi\in[0,2\pi].$ Since $\mathcal{L}$ is a smooth, egg-shaped contour, the solution is unique. Clearly, $\rho(\phi)=\frac{\delta(\phi)}{\cos(\phi)}$, and the parametrization of $\mathcal{L}$ is fully specified. Then, the mapping from $z\in\mathcal{C}_{z}^{+}$ to $y\in\mathcal{L}^{+}$, where $z=e^{i\phi}$ and $y=\rho(\psi(\phi))e^{i\psi(\phi)}$, satisfying $\gamma_{0}(0)=0$, $\gamma_{0}(z)=\overline{\gamma_{0}(z)}$ is uniquely determined by,

[TABLE]

i.e., $\psi(.)$ is uniquely determined as the solution of a Theodorsen integral equation with $\psi(\phi)=2\pi-\psi(2\pi-\phi)$. Due to the correspondence-boundaries theorem, $\gamma_{0}(z)$ is continuous in $\mathcal{C}_{z}\cup\mathcal{C}_{z}^{+}$.

Appendix 0.C On the analyticity of $\Pi_{j}(x,y)$ close to $p=0$

We focus only on the analyticity of $\Pi_{0}(x,y)$ in a neighborhood of $p=0$ by using a variant of the implicit function theorem on the functional equation (6). The analyticity of $\Pi_{1}(x,y)$, follows directly by the analyticity of $\Pi_{0}(x,y)$ from the second in (4). We follow the lines in [39, 15], and use the implicit function theorem for Banach spaces (see Theorem 10.2.3, p. 272 in [12]). Define the mapping $f:S\subset\mathbb{C}\times B_{2}\to B_{3}\times\mathbb{C}$,

[TABLE]

where $S$ contains the point $(0,V_{0}^{(0)})$, $H$, $A$, $B$, $C$, are as in (20) and (5) respectively, $B_{2}$ be the Banach space comprising all bivariate analytic bounded functions in $\mathbb{D}^{2}$, with $\mathbb{D}$ the open complex unit disk, and $B_{3}$ be the Banach space comprising all trivariate analytic bounded functions in $\mathbb{D}^{3}$ that have a limit of 0 for the first two arguments going to 1.

Since $H$, $A$, $B$, $C$ are bounded analytic functions in $\mathbb{D}^{3}$, and since $f$ is affine in $\Pi_{0}$ and $p$, it is easily seen that $f$ is $r$-times continuously differentiable for all $r$. Note also that $f(0,V_{0}^{(0)})=[0,0]$. Then, the (Banach space) derivative of $f$ at the point $(0,V_{0}^{(0)})$ [12] equals

[TABLE]

We need to show that this mapping is a homeomorphism. Indeed,

$df(0,V_{0}^{(0)})$ is a continuous mapping for the same reasons that the mapping $f$ itself is continuous. 2. 2.

For given $\Pi_{0}^{(1)}$, $\Pi^{(2)}_{0}$, let $df(0,V_{0}^{(0)})(\Pi_{0}^{(1)})=df(0,V_{0}^{(0)})(\Pi_{0}^{(1)})$. Then,

[TABLE]

or equivalently $f(0,\Pi_{0}^{(1)}-\Pi_{0}^{(2)})=(0,-\frac{\tau}{\tau+\gamma})$, which in turn has the zero solution as a unique solution [3], and thus $\Pi^{(0)}_{1}=\Pi^{(0)}_{2}$ so that $df(0,V_{0}^{(0)})$ is injective. 3. 3.

To show that $df(0,V_{0}^{(0)})$ is surjective, we solve the $df(0,V_{0}^{(0)})(\Pi_{0})=(g,c)$ with $g$ a bivariate analytic bounded function in $\mathbb{D}^{2}$ with limit 0 for its arguments going to 1, and $c$ a complex number. The solution is

[TABLE] 4. 4.

The $\Pi_{0}$ obtained previously equals $(df(0,V_{0}^{(0)}))^{-1}$, which is readily seen that it is continuous.

Thus, $\Pi_{0}\to df(0,V_{0}^{(0)})(\Pi^{(0)})$ is a linear homeomorphism and using Theorem 10.2.3 in [12], $\Pi_{0}(x,y)$ is $r$-times differentiable at $p=0$. Having this result, and using the second in (4), $\Pi_{1}(x,y)$ is also $r$-times differentiable at $p=0$.