A Two-Server Queue with Interdependence between Arrival and Service Processes

Abstract

In this paper, we analyse a queueing system with two servers where the arrival and service processes are interdependent. The evolution of these processes is governed by transitions on the product space of three Markov chains, which are descriptors of the arrival and service processes. The transitions in this Markov chain follow a semi-Markov rule and the exponential distribution governs the sojourn times in the states. The stability condition of the system is derived and the stationary distribution is calculated for the system in equilibrium. Several important performance measures are provided, and numerical illustrations of the model are presented.

1. Introduction

Queueing theory is a branch of applied mathematics that studies waiting lines or queues. It finds widespread applications in various fields, including telecommunications, computer networks, manufacturing systems, and service industries. Two server queueing models are particularly significant, as they offer insights into scenarios where multiple servers handle arriving customers or tasks. While traditional queueing models assume independent arrival and service processes, in many practical situations these processes are interconnected in a way that can significantly impact system performance and behaviour. The study of queues with interdependent arrival and service processes helps in understanding how these mutual influences shape queue dynamics, waiting times, resource utilisation, and overall system efficiency.

1.1. Literature Review of Queues with Interdependence in Arrival and Service Processes

Mitchell et al. [1] investigated an M/M/1 queue in which the arrival and service processes are linked through correlated random variables, conforming to a bivariate exponential distribution. This approach sheds light on how correlations between customer service times and the time between their arrival impact the system. Courtois et al. [2] generalized the M/G/1 queuing process by assuming both the arrival and service rates to be essentially arbitrary functions of the number of customers currently using the system. In addition, they assumed that the amount of service demanded by a customer is conditioned by the queue length at the moment service is begun for that customer. Sengupta [3] investigated a single-server semi-Markovian queue with a first-come-first-served policy in which both the arrival and service processes follow a semi-Markov process. Boxma et al. [4] considered a storage model with a dual interpretation: as a distinctive queueing model featuring interdependence between service requests and subsequent inter-arrival times, and as a fluid production/inventory model operating within a two-state random environment.

Muller et al. [5] studied a queueing system characterised by partial correlation, focusing on the interplay between the quantity of work (service time) introduced by a customer and the subsequent inter-arrival time, which the authors assumed to be interdependent. The focus of the study revolved around a distinctive single-server multi-station alternating queue configuration in which the interplay between preparation and service times involves intricate auto- and cross-correlation patterns. Their investigation was carried out across two distinct scenarios. Vlasiou et al. [6] addressed a configuration involving a single-server multi-station alternating queue characterised by interrelated preparation and service times and exhibiting both auto- and cross-correlations. Adan et al. [7] presented a comprehensive investigation into a single-server queue system in which both the inter-arrival times and the service times depend on a common discrete time Markov chain. This model represents an extension of the well-known MAP/G/1 queue, incorporating interdependencies between inter-arrival and service times. Iyer et al. [8] analysed queueing models characterised by joint density of the inter-arrival and service times. These models involve a mixture of joint densities, and find natural applications in scenarios involving a single server catering to multiclass populations through a shared queue. Additionally, this approach can help to model interdependencies between inter-arrival and service times in the context of queue control models. Buchholz et al. [9] focused on a queueing system in which the inter-arrival times exhibit correlation and the service times are correlated with the inter-arrival times. They established a compelling connection between this correlated queueing model and the MMAP[K]/PH[K]/1 queue, a well-studied structure to which matrix geometric algorithms are readily applicable. Kim et al. [10] analysed the wait time distribution in an M/M/1 queue in which the inter-arrival time between the n-th and $(n+1)$-th customers along with the service time of the n-th customer are modelled as correlated random variables characterised by Downton’s bivariate exponential distribution. Krishnamoorthy et al. [11] examined a queueing system with finite Markov-dependent arrival and service batch sizes, and additionally considered correlated successive inter-arrival and service time durations. Dai et al. [12] investigated a unique category of correlated queue in which the service duration of an arriving customer hinges upon the inter-arrival time separating them from the previous customer.

Moiseev et al. [13] examined a two-stage queuing system with infinite servers in which customers arrive according to a Markovian Arrival Process (MAP). In this system, the service times at the first stage and the probability of transitioning to the second stage depend on the type of request, which is determined by the state of the underlying Markov chain of MAP when the request arrives. The authors’ analysis of this model utilized the multi-dimensional dynamic screening method.

1.2. Literature Review of Queues with Multiservers

The literature abounds with extensive research concerning queueing models that involve multiple servers and operate under the assumption of independent arrival and service processes. Kleinrock [14] laid the foundation for analysing queueing systems with multiple servers. The earliest works focused on fundamental models characterised by independent customer arrivals and exponential service durations. These pioneering studies offered essential insights into the behaviour of queueing systems, particularly in the context of systems featuring multiple servers. These foundational investigations laid the groundwork for subsequent advancements in the field. Cohen [15] played a crucial role in advancing this trajectory by expanding the analytical framework to accommodate batch arrivals and service times. This extension was a pivotal step, contributing to a more comprehensive understanding of the dynamics of queueing systems with multiple servers. Neuts et al. [16] addressed a notable challenge by highlighting the analytical complexity of queueing systems involving more than two heterogeneous servers. They aptly observed that in such scenarios the steady-state behaviour of the system can only be explored through algorithmic approaches. Building upon this observation, Krishna Kumar et al. [17] studied an M/M/2 queueing system characterised by heterogeneous servers.

There are a number of recent papers dealing with two-server queueing systems, including [18,19,20,21,22,23,24,25,26,27,28], to which interested readers may refer.

A semi-Markov process is both a Markov renewal process and a generalization of the Markov process. Unlike Markov renewal processes, which focus on the number of visits to each state, a semi-Markov process considers random variables describing the time spent in each state before transitioning. In discrete-time Markov processes, transitions occur after fixed unit time intervals; in continuous-time Markov processes, transitions follow exponential distributions for the time spent in each state. A semi-Markov process resembles a Markov process in terms of state-to-state transitions; however, it is set apart by the time spent in each state before moving to the next state being a random variable with parameter(s) that depend on both the state it is currently in and the one to be visited next.

In a semi-Markov process, changes of the occupied states take place according to a Markov chain rule; nevertheless, a random amount of time between changes is required. Suppose that a process can be in any one of N states $1,2,3\dots ,N$ and that each time it enters state $i,i=1,2,3,\dots ,N$ it remains there for a random amount of time before transitioning to state j with probability $P_{ij}$. The state j that is to be visited next is chosen at the time that state i is visited. The sojourn time in i depends on both i and j. Let the time until the transition from i to j have a distribution $F_{ij}$, and let $Z\left(t\right)$ denote the state at time t; then, $\left\{Z\right(t):t\ge 0\}$ is called a semi-Markov processes.

To define a semi-Markov process, it is necessary to have a finite state space Markov chain (MC). This MC can be one-dimensional or higher, and must be finite. Associated with the MC is its initial probability vector (IPV) and one-step transition probability matrix (TPM). The semi-Markov process evolves as follows. First, choose an initial state j according to the IPV. A particle occupies this state. At this moment, the next state (say, k) to be visited by the particle is chosen; this is sampled using the one-step TPM. The sojourn time of the particle in j is distributed based on the mean time in j as ${\mu }_{jk}$. At the end of this time, the particle transits to k; at this epoch, the next state (say, m) to be visited is chosen using the one-step TPM. In k, the sojourn time has a mean value ${\mu }_{km}$. This process continues; the continuous time process is referred to as a semi-Markov process. The sojourn time distribution $F_{ij}$ can be exponential or any general distribution (with the non-negative part of the real line as support). In this paper, the sojourn time distribution is taken as the exponential distribution.

We consider a two-server queueing system in which the arrival and service processes are interdependent. Numerous studies in the existing literature have investigated the interdependence among random variables and processes through various analytical methods. In this paper, we analyse the interdependence among arrival and service processes using a semi-Markov approach. This is a new approach first introduced by Achyutha Krishnamoorthy in 2021 (see [29,30]). The procedure is as follows. Suppose that n processes evolve in an interdependent fashion, as explained below. Assume that these n processes are generated through transitions in the states of n finite state space-distinct Markov chains. They may or may not have absorbing states. Considering the product space of these Markov chains, a Markov chain is imposed on this product space. This chain has an initial probability vector and a one-step transition probability matrix. The combined process evolves through transitions in the elements of the product space. We may assume that these transitions are governed by a semi-Markov rule, with the exception of the sojourn time in each stage|state (n-tuples) being exponentially distributed based on a parameter with dependence on the current state and the one to be visited next (as governed by the one-step transition probability matrix on the product space Markov chain).

In [31], the authors examined two distinct queueing models. Model I used a single-server working vacation queueing system in which the arrival and service processes were interdependent. The evolution of the arrival and service processes occurs through transitions within the product space of two separate Markov chains of orders m and n, respectively. The transitions in the product spaces of these Markov chains are governed by a semi-Markov rule, with sojourn times in the states governed by an exponential distribution with parameters depending on the state it is in and the state to be visited next, which are determined according to the Markov chain rule. During working vacation periods, service is delivered at a reduced rate. The duration of these vacation intervals follows an exponential distribution with the parameter $\eta$. In contrast, the second model is based on independent arrival and service processes following phase-type distributions with representations ($\mathit{\alpha }$, T) of order m and ($\mathit{\beta }$, $S)$ of order n, respectively. The service time during normal operation is the phase-type distribution indicated above, whereas that during working vacation is a phase-type distribution with representation ($\mathit{\beta }$, $\theta S),0<\theta <1$.

In the present paper, we analyse a queueing system with two servers having interdependent arrival and service processes. The evolution of these processes is governed by transitions in the product space of three Markov chains. The transitions in this Markov chain follow a semi-Markov rule and the exponential distribution governs the sojourn times in the states.

1.3. Motivation

In queueing theory literature, many papers have been developed to analyse arrival and service processes which are mutually independent with the exception that service can be rendered only when the server and at least one customer to receive the service are available. A few papers have considered successive inter-arrival times and/or successive service times to be correlated, for example, Markovian or batch Markovian arrival/ service processes. Even fewer papers have considered the interdependence of arrival and service processes. For example, the inter-arrival times $A_n$ between the $n-1$ and n-th customers and the service time $S_n$ of the n-th customer were correlated for $n=1,2,\dots .$ by Sengupta [3,32]. Van Houdt [33] subsequently filled certain gaps in the work of Sengupta [3]. In all these papers, the authors made specific distributional assumptions with respect to the arrival and service times. Thus, the question arises as to whether it is necessary to specify the distributions of the inter-arrival and service times of customers. It is this question that led us to design the presented approach to multi-server queues.

Thus, in this paper we make no assumptions about the inter-arrival and service time distribution. More importantly, no such assumption can be meaningfully made in the context of modelling this interdependence. There are at most three papers [29,30,31] that have appeared in the stochastic modelling literature employing this new approach, and no cases involving more than one server have been considered under this framework in the extant literature. Thus, the methodology presented here is quite novel. This procedure can be extended to queues of a finite number with more than two servers. The dimension increases by one as the number of servers increases.

Markov chains provide the simplest dependence among a sequence of random variables. In the present modelling, we assume that arrival and service rendered by the two servers evolve within the product space of three Markov chains: one each for arrival, service by server 1, and service by server 2. Thus, our approach is completely new, especially in service systems with multiple servers.

1.4. Practical Applications

In the correlated inter-arrival time and service time, we can use the following example: with $A_n$ being the sequence of independent and identically distributed inter-arrival times and $S_n$ the sequence of service times of successive customers, the stress is on the correlation, that is, whether the correlation is positive, negative, or zero. A negative correlation indicates that if the inter-arrival time $A_n$ is large, the corresponding service time $S_n$ turns out to be comparatively small, and vice versa. In the positive correlation case, this indicates a “direct proportion” between the two. In our model of interdependence, quick transitions in the arrival phases indicate faster arrival. Accordingly, the servers can accelerate their service speed or, equivalently, their service rate. In the correlated arrival and service processes indicated above, the n-th customer to arrive receives service at a future time after arrival, as the server can be busy at the customer’s arrival epoch. While the situation in our model is the same, the reflection of the interdependence of arrival and service is seen simultaneously. This helps to accelerate or decelerate the rate of service provided by the servers. These kinds of adjustments can be seen in the production process; if customers arrive at a fast rate, the production needs to be adjusted to ensure that the queue size does not blow up.

In large airports, when non-domestic flights flights land the passengers are required to proceed to the immigration section. A queue of such passengers is formed; at the head of the queue, the person who manages the distribution of passengers directs them to different immigration counters depending on the slots available at each counter. The processing of each passenger passes through distinct stages in a forward direction. With only slight modifications, the model discussed in the paper can be implemented to clear long queues of passengers efficiently.

Another practical application of a two-server queueing model with interdependent arrival and service processes can be found in scenarios involving distributed data processing or computing. In such scenarios, data packets or tasks arrive at a network of servers, and the rate of arrival and processing time can be influenced by the current state of the servers in the network. The arrival rate of data packets may be influenced by the current number of idle or busy servers. For example, if both servers are idle, the arrival rate may increase, as more data processing capacity is available. Conversely, if both servers are busy, the arrival rate might decrease as new data packets wait to be processed. In such a scenario, the service rate for each server can depend on its current workload.

The salient features of the model discussed in this paper are as follows:

• We introduce a new approach to analysing two-server queueing systems with interdependent arrival and service processes.
• The presented analysis does not involve the Kronecker product/sum, unlike the queueing models with independent arrival and service processes analysed using matrix geometric methods.

Notations and abbreviations used in this paper:

• $CTMC$: continuous-time Markov chain.
• $LIQBD$: level-independent quasi-birth and death.
• $\mathit{e}$: column vector of 1s of appropriate order.
• $QBD$: quasi-birth and death.

The remainder of the paper is organized as follows: in Section 2, the model is mathematically formulated; in Section 3, we perform a steady-state analysis of the studied queuing model after establishing the stability conditions of the system; in Section 4, performance measures are computed and presented; finally, numerical illustrations of the two models are discussed in Section 5.

2. Model Description and Mathematical Formulation

Consider a queueing system with two heterogeneous servers $S_1$ and $S_2$ in which the arrival and service processes are interdependent. When a customer arrives with both servers idle, server $S_1$ provides service. When both servers are occupied, an arriving customer joins an infinite capacity queue. When the system becomes empty, both servers remain idle. The arrival process passes through several stages, including back and forth. An arriving customer who finds a free server enters for service immediately. The service process provided by server $S_1$ has $n_1$ stages, while that provided by server $S_2$ has $n_2$ stages. These proceed in the forward direction. The arrival process has m stages. Three Markov chains govern the arrival and service processes. Service provided by $S_1$ is governed by the Markov chain $\mathcal{X}=\left\{X_n\right\}$ and service provided by $S_2$ is governed by the Markov chain $\mathcal{Y}=\left\{Y_n\right\}$. The state space of the Markov chain $\left\{X_n\right\}$ is $\{1,2,3,\dots ,n_1,n_1+1\}$, while that of the Markov chain $\left\{Y_n\right\}$ is $\{1,2,3,\dots ,n_2,n_2+1\}$. The arrival process is governed by a Markov chain $\mathcal{Z}=\left\{Z_n\right\}$ with state space $\{1,2,3,\dots ,m,m+1\}$. We can consider the Markov chain on the product space $\mathcal{S}=\mathcal{X}\times \mathcal{Y}\times \mathcal{Z}$ with state space $\{\{1,2,3,\dots ,n_1,n_1+1\}\times$$\{1,2,3,\dots ,n_2,n_2+1\}\times \{1,2,3,\dots ,m,m+1\}\}$.

The absorbing states of this Markov chain are $\{(i,j,m+1):1\le i\le n_1;1\le j\le n_2\}\bigcup$ $$\begin{gathered} \{(i,n_2+1,k):1\le i\le n_1;1\le k\le m\}\bigcup \{(n_1+1,j,k):1\le j\le n_2;1\le k\le m\}\bigcup \{(n_1+1, \\ {n}_2+1,k):1\le k\le m\}\bigcup \{(n_1+1,j,m+1):1\le j\le n_2\}\bigcup \end{gathered}$$ $\{(i,n_2+1,m+1):1\le i\le n_1\}$ $\bigcup \left\{(n_1+1,n_2+1,m+1)\right\}$. Changes in the first coordinate due to transitions indicate service phase changes provided by server $S_1$, those in the second coordinate indicate service phase changes provided by server $S_2$, and those in the third coordinate indicate arrival phase changes. The transitions are interdependent in the sense that the sojourn time in any state $(i,j,k)$ depends on this state as well as the one (say, $(i^{{}^{\prime }},j^{{}^{\prime }},k^{{}^{\prime }})$) to be visited next. This sojourn time is assumed to be exponentially distributed with parameter ${\delta }_{(i,j,k)(i^{{}^{\prime }},j^{{}^{\prime }},k^{{}^{\prime }})}$. Because the transitions are interdependent, within a short interval it is possible for none, one, two, or even all coordinates to change with positive probability when a transition occurs. The state space indicated above is for this general case; however, we assume here that at most one change takes place with positive probability. This assumption leads to an infinitesimal generator which is highly sparse. Let the initial probability vector of the arrival process be $\mathit{\alpha }$$=($$\overline{\mathit{\alpha }}$$,{\alpha }_{m+1})$, where $\overline{\mathit{\alpha }}$$=({\alpha }_1,{\alpha }_2,{\alpha }_3,\dots {\alpha }_m)$. We can further simplify the assumption that the service processes for both servers start from stage 1 and move in the forward direction only, that is, the service processes are in the order $1\to 2\to 3\to 4\cdots n_1\to n_1+1$ and $1\to 2\to 3\to 4\cdots n_2\to n_2+1$, respectively. The initial probability vectors of the service processes of both the servers are $(1,0,0,\dots )$$n_1+1$ and $n_2+1$ component vectors, respectively. Further, to ensure that arrivals do not occur too quickly, we assume that ${\alpha }_{m+1}=0$.

2.1. The QBD Process

The model described in Section 2 can be studied as a Level-Independent Quasi-Birth–Death (LIQBD) process (see Latouche et al. [34]). First, we define the following notations. At time t, let:

• $N\left(t\right)$ be the number of customers in the system.
• $S_1\left(t\right)$ be the the phase of the service provided by server $S_1$.
• $S_2\left(t\right)$ be the phase of the service provided by server $S_2$.
• $A\left(t\right)$ be the phase of arrival.

Here, $\{(N\left(t\right),S_1\left(t\right),S_2\left(t\right),A\left(t\right)):t\ge 0\}$ is an LIQBD process with state space

$\tilde{\Omega }$= $$\begin{gathered} \{\{(0,\ast ,\ast ,k):1\le k\le m\}\bigcup \{(1,i,\ast ,k):1\le i\le n_1;1\le k\le m\}\bigcup \{(1,\ast ,j,k):1\le j\le n_2;1\le k\le m\} \\ \bigcup \{(q,i,j,k):1\le i\le n_1;1\le j\le n_2;1\le k\le m;q\ge 2\}\} \end{gathered}$$. In the absence of customers, no service can be provided; this is indicated by the * in the position of the service coordinates (the third and fourth coordinates of the 4-tuple). If only one customer is in the system, either $S_1$ or $S_2$ is idle, when $S_2$ is idle it is represented by $(1,i,\ast ,k)$ and when $S_1$ is idle it is represented by $(1,\ast ,j,k)$, where $1\le i\le n_1;1\le j\le n_2;1\le k\le m$.

2.2. Transitions

The transitions are described in

Table 1. Transition rates.

From	To	Rate		Remarks
$(0,\ast ,\ast ,k)$	$(0,\ast ,\ast ,k^{{}^{\prime }})$	${\delta }_{(\ast ,\ast ,k)(\ast ,\ast ,k^{{}^{\prime }})}$	$1\le k,k^{{}^{\prime }}\le m,$	arrival phase change
$(0,\ast ,\ast ,k)$	$(1,1,\ast ,k^{{}^{\prime }})$	${\alpha }_{k^{{}^{\prime }}}{\delta }_{(\ast ,\ast ,k)(1,\ast ,m+1)}$	$1\le k,k^{{}^{\prime }}\le m$	arrival occurs
$(1,i,\ast ,k)$	$(1,i,\ast ,k^{{}^{\prime }})$	${\delta }_{(i,\ast ,k)(i,\ast ,k^{{}^{\prime }})}$	$1\le k,k^{{}^{\prime }}\le m;1\le i\le n_1$	arrival phase change
$(1,n_1,\ast ,k)$	$(0,\ast ,\ast ,k)$	${\delta }_{(n_1,\ast ,k)(n_1+1,\ast ,k)}$	$1\le k\le m$	service completion at $S_1$
$(1,i,\ast ,k)$	$(1,i+1,\ast ,k)$	${\delta }_{(i,\ast ,k)(i+1,\ast ,k)}$	$1\le k\le m;1\le i\le n_1-1$	service phase change at $S_1$
$(1,i,\ast ,k)$	$(2,i,1,k^{{}^{\prime }})$	${\delta }_{(i,\ast ,k)(i,1,m+1)}{\alpha }_{k^{{}^{\prime }}}$	$1\le k,k^{{}^{\prime }}\le m,1\le i\le n_1$	arrival occurs
$(1,\ast ,j,k)$	$(2,1,j,k^{{}^{\prime }})$	${\delta }_{(\ast ,j,k)(1,j,m+1)}{\alpha }_{k^{{}^{\prime }}}$	$1\le k,k^{{}^{\prime }}\le m,1\le j\le n_2$	arrival occurs
$(1,\ast ,j,k)$	$(1,\ast ,j+1,k)$	${\delta }_{(\ast ,j,k)(\ast ,j+1,k)}$	$1\le j\le n_1-1$; $1\le k\le m;$	service phase change at $S_2$
$(1,\ast ,j,k)$	$(1,\ast ,j,k^{{}^{\prime }})$	${\delta }_{(\ast ,j,k)(\ast ,j,k^{{}^{\prime }})}$	$1\le j\le n_2;1\le k,k^{{}^{\prime }}\le m$	arrival phase change
$(1,\ast ,n_2,k)$	$(0,\ast ,\ast ,k)$	${\delta }_{(\ast ,n_2,k)(\ast ,n_2+1,k)}$	$1\le k\le m$	service completion at $S_2$
$(2,n_1,j,k)$	$(1,\ast ,j,k)$	${\delta }_{(n_1,j,k)(n_1+1,j,k)}$	$1\le j\le n_2;1\le k\le m$	service completion at $S_1$
$(2,i,n_2,k)$	$(1,i,\ast ,k)$	${\delta }_{(i,n_2,k)(i,n_2+1,k)}$	$1\le i\le n_1;1\le k\le m$	service completion at $S_2$
$(q,i,j,k)$	$(q+1,i,j,k^{{}^{\prime }})$	${\alpha }_{k^{{}^{\prime }}}{\delta }_{(i,j,k)(i,j,m+1)}$	$1\le k,k^{{}^{\prime }}\le m;1\le i\le n_1;1\le j\le n_2;q\ge 2$	arrival occurs
$(q,i,j,k)$	$(q,i,j,k^{{}^{\prime }})$	${\delta }_{(i,j,k)(i,j,k^{{}^{\prime }})}$	$1\le k,k^{{}^{\prime }}\le m;1\le i\le n_1,1\le j\le n_2;q\ge 2$	arrival phase change
$(q,i,j,k)$	$(q,i+1,j,k)$	${\delta }_{(i,j,k)(i+1,j,k)}$	$1\le k\le m;1\le i\le n_1-1;1\le j\le n_2;q\ge 2$	service phase change at $S_1$
$(q,i,j,k)$	$(q,i,j+1,k)$	${\delta }_{(i,j,k)(i,j+1,k)}$	$1\le k\le m;1\le i\le n_1;1\le j\le n_2-1;q\ge 2$	service phase change at $S_2$
$(q,n_1,j,k)$	$(q-1,1,j,k)$	${\delta }_{(n_1,j,k)(n_1+1,j,k)}$	$1\le k\le m;1\le j\le n_2;q\ge 3$	service completion at $S_1$
$(q,i,n_2,k)$	$(q-1,i,1,k)$	${\delta }_{(i,n_2,k)(i,n_2+1,k)}$	$1\le k\le m;1\le i\le n_1;q\ge 3$	service completion at $S_2$

.

The infinitesimal generator of the CTMC is

$$Q^{\ast }=\left[\begin{array}{cccc}{\mathcal{B}}_1& {\mathcal{B}}_0& & \\ {\mathcal{B}}_2& {\mathcal{C}}_1& {\mathcal{C}}_0& \\ & {\mathcal{C}}_2& {\mathcal{A}}_1& {\mathcal{A}}_0& \\ & & {\mathcal{A}}_2& {\mathcal{A}}_1& {\mathcal{A}}_0& \\ & & & {\mathcal{A}}_2& {\mathcal{A}}_1& {\mathcal{A}}_0& \\ & & & \ddots & \ddots & \ddots \end{array}\right].$$

Here, ${\mathcal{B}}_1$ is a square matrix of order m which contains the transition rates within level 0, ${\mathcal{B}}_0$ is a matrix of order $m\times (n_1+n_2)m$ which contains transition rates from level 0 to level 1, ${\mathcal{B}}_2$ is a matrix of order $(n_1+n_2)m\times m$ which contains transition rates from level 1 to level 0, ${\mathcal{C}}_1$ is a square matrix of order $(n_1+n_2)m$ which contains the transition rates within the level 1, ${\mathcal{C}}_0$ is a matrix of order $(n_1+n_2)m\times n_1{n}_{2}m$ which contains transition rates from level 1 to level 2, ${\mathcal{C}}_2$ is a matrix of order $n_1{n}_{2}m\times (n_1+n_2)m$ which contains transition rates from level 2 to level 1, ${\mathcal{A}}_0$ represents transition rates from level n to level $n+1$ for $n\ge 2$, ${\mathcal{A}}_1$ represents transition rates within the level n for $n\ge 2$, and ${\mathcal{A}}_2$ represents transition rates from the level n to level $n-1$ for $n\ge 3$. All of these are square matrices of order $n_1{n}_{2}m$.

The matrix ${\mathcal{B}}_1$ is provided by

$$\left[\begin{array}{ccccc}{b}_1& {\delta }_{(\ast ,\ast ,1)(\ast ,\ast ,2)}& {\delta }_{(\ast ,\ast ,1)(\ast ,\ast ,3)}& \cdots \cdots & {\delta }_{(\ast ,\ast ,1)(\ast ,\ast ,m)}\\ {\delta }_{(\ast ,\ast ,2)(\ast ,\ast ,1)}& b_2& {\delta }_{(\ast ,\ast ,2)(\ast ,\ast ,3)}& \cdots \cdots & {\delta }_{(\ast ,\ast ,2)(\ast ,\ast ,m)}\\ {\delta }_{(\ast ,\ast ,3)(\ast ,\ast ,1)}& {\delta }_{(\ast ,\ast ,3)(\ast ,\ast ,2)}& b_3& \cdots \cdots & {\delta }_{(\ast ,\ast ,3)(\ast ,\ast ,m)}\\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ {\delta }_{(\ast ,\ast ,m)(\ast ,\ast ,1)}& {\delta }_{(\ast ,\ast ,m)(\ast ,\ast ,2)}& {\delta }_{(\ast ,\ast ,m)(\ast ,\ast ,3)}& \cdots \cdots & b_m\end{array}\right],$$

which is a square matrix of order m where $b_i$ = $-[{\displaystyle \sum _{j=1,j\ne i}^m}{\delta }_{(\ast ,\ast ,i)(\ast ,\ast ,j)}+{\delta }_{(\ast ,\ast ,i)(1,\ast ,m+1)}]$; $1\le i\le m$.

We define the matrix F as

$$\left[\begin{array}{ccccc}{\alpha }_1{\delta }_{(\ast ,\ast ,1)(1,\ast ,m+1)}& {\alpha }_2{\delta }_{(\ast ,\ast ,1)(1,\ast ,m+1)}& {\alpha }_3{\delta }_{(\ast ,\ast ,,1)(1,\ast ,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(\ast ,\ast ,1)(1,\ast ,m+1)}\\ {\alpha }_1{\delta }_{(\ast ,\ast ,2)(1,\ast ,m+1)}& {\alpha }_2{\delta }_{(\ast ,\ast ,2)(1,\ast ,m+1)}& {\alpha }_3{\delta }_{(\ast ,\ast ,2)(1,\ast ,m+)}& \cdots \cdots & {\alpha }_m{\delta }_{(\ast ,\ast ,2)(1,\ast ,m+1)}\\ {\alpha }_1{\delta }_{(\ast ,\ast ,3)(1,\ast ,m+1)}& {\alpha }_2{\delta }_{(\ast ,\ast ,3)(1,\ast ,m+1)}& {\alpha }_3{\delta }_{(\ast ,\ast ,3)(1,\ast ,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(\ast ,\ast ,3)(1,\ast ,m+1)}\\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ {\alpha }_1{\delta }_{(\ast ,\ast ,m)(1,\ast ,m+1)}& {\alpha }_2{\delta }_{(\ast ,\ast ,m)(1,\ast ,m+1)}& {\alpha }_3{\delta }_{(\ast ,\ast ,m)(1,\ast ,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(\ast ,\ast ,m)(1,\ast ,m+1)}\end{array}\right],$$

which is a square matrix of order m.

Then, ${\mathcal{B}}_0$ can be expressed as $\left[\begin{array}{ccccc}F& \mathbf{0}& \mathbf{0}& \cdots \cdots & \mathbf{0}\end{array}\right]$.

We write $D_1$ as

$$\left[\begin{array}{cccc}{\delta }_{(n_1,\ast ,1)(n_1+1,\ast ,1)}& & & \\ & {\delta }_{(n_1,\ast ,2)(n_1+1,\ast ,2)}& & & \\ & & {\delta }_{(n_1,\ast ,3)(n_1+1,\ast ,3)}& & \\ & & & \ddots \\ & & & & {\delta }_{(n_1,\ast ,m)(n_1+1,\ast ,m)}\end{array}\right],$$

which is a square matrix of order m.

We define the matrix $D_2$ as

$$\left[\begin{array}{cccc}{\delta }_{(\ast ,n_2,1)(\ast ,n_2+1,1)}& & & \\ & {\delta }_{(\ast ,n_2,2)(\ast ,n_2+1,2)}& & & \\ & & {\delta }_{(\ast ,n_2,3)(\ast ,n_2+1,3)}& & \\ & & & \ddots \\ & & & & {\delta }_{(\ast ,n_2,m)(\ast ,n_2+1,m)}\end{array}\right],$$

which is a square matrix of order m.

Then, ${\mathcal{B}}_2=\left[\begin{array}{c}\mathbf{0}\\ \mathbf{0}\\ ⋮\\ D_1\\ \mathbf{0}\\ \mathbf{0}\\ ⋮\\ D_2\end{array}\right]$.

We define

$$E_i\left[\begin{array}{ccccc}{e}_{i1}& {\delta }_{(i,\ast ,1)(i,\ast ,2)}& {\delta }_{(i,\ast ,1)(i,\ast ,3)}& \cdots & {\delta }_{(i,\ast ,1)(i,\ast ,m)}\\ {\delta }_{(i,\ast ,2)(i,\ast ,1)}& e_{i2}& {\delta }_{(i,\ast ,2)(i,\ast ,3)}& \cdots & {\delta }_{(i,\ast ,2)(i,\ast ,m)}\\ {\delta }_{(i,\ast ,3)(i,\ast ,1)}& {\delta }_{(i,\ast ,3)(i,\ast ,2)}& e_{i3}& \cdots & {\delta }_{(i,\ast ,3)(i,\ast ,m)}\\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ {\delta }_{(i,\ast ,m)(i,\ast ,1)}& {\delta }_{(i,\ast ,m)(i,\ast ,2)}& {\delta }_{(i,\ast ,m)(i,\ast ,3)}& \cdots & e_{im}\end{array}\right]$$

which as a square matrix of order m, where

$$\begin{array}{c}\hfill e_{iq}=-[\sum _{k=1,q\ne k}^m{\delta }_{(i,\ast ,q)(i,\ast ,k)}+{\delta }_{(i,\ast ,q)(i,1,m+1)}+{\delta }_{(i,\ast ,q)(i+1,\ast ,q)}];1\le i\le n_1,1\le q\le m;\end{array}$$

$$G_i=\left[\begin{array}{cccc}{\delta }_{(i,\ast ,1)(i+1,\ast ,1)}& & & \\ & {\delta }_{(i,\ast ,2)(i+1,\ast ,2)}& & & \\ & & {\delta }_{(i,\ast ,3)(i+1,\ast ,3)}& & \\ & & & \ddots \\ & & & & {\delta }_{(i,\ast ,m)(i+1,\ast ,m)}\end{array}\right]$$

as a square matrix of order m with $1\le i\le n_1-1$;

$$H_j=\left[\begin{array}{ccccc}{h}_{j1}& {\delta }_{(\ast ,j,1)(\ast ,j,2)}& {\delta }_{(\ast ,j,1)(\ast ,j,3)}& \cdots & {\delta }_{(\ast ,j,1)(\ast ,j,m)}\\ {\delta }_{(\ast ,j,2)(\ast ,j,1)}& h_{j2}& {\delta }_{(\ast ,j,2)(\ast ,j,3)}& \cdots & {\delta }_{(\ast ,j,2)(\ast ,j,m)}\\ {\delta }_{(\ast ,j,3)(\ast ,j,1)}& {\delta }_{(\ast ,j,3)(\ast ,j,2)}& h_{j3}& \cdots & {\delta }_{(\ast ,j,3)(\ast ,j,m)}\\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ {\delta }_{(\ast ,j,m)(\ast ,j,1)}& {\delta }_{(\ast ,j,m)(\ast ,j,2)}& {\delta }_{(\ast ,j,m)(\ast ,j,3)}& \cdots & h_{jm}\end{array}\right]$$

as a square matrix of order m, where

$$\begin{array}{c}\hfill h_{jp}=-[\sum _{p=1,p\ne k}^m{\delta }_{(\ast ,j,p)(\ast ,j,k)}+{\delta }_{(\ast ,j,p)(\ast ,j,m+1)}+{\delta }_{(\ast ,j,p)(\ast ,j+1,p)}];1\le j\le n_2,1\le k\le m,\end{array}$$

and

$$J_j=\left[\begin{array}{cccc}{\delta }_{(\ast ,i,1)(\ast ,i+1,1)}& & & \\ & {\delta }_{(\ast ,i,2)(\ast ,i+1,2)}& & & \\ & & {\delta }_{(\ast ,i,3)(\ast ,i+1,3)}& & \\ & & & \ddots \\ & & & & {\delta }_{(\ast ,i,m)(\ast ,i+1,m)}\end{array}\right]$$

as a square matrix of order m with $1\le j\le n_2-1$.

Using these, we can write

$${\mathcal{C}}_1=\left[\begin{array}{ccccccc}{E}_1& G_1& & & & & \\ & E_2& G_2& & & & \\ & & \ddots \ddots \\ & & & E_{n_1-1}& G_{n_1-1}& & & \\ & & & & E_{n_1}& \mathbf{0}& \mathbf{0}\\ & & & & \mathbf{0}& H_1& J_1\\ & & & & & & H_2& J_2\\ & & & & & & & \ddots \ddots & \\ & & & & & & & & & H_{n_2-1}& J_{n_2-1}& \\ & & & & & & & & & & H_{n_2}\end{array}\right].$$

Let

$$K_i=\left[\begin{array}{ccccc}{\alpha }_1{\delta }_{(i,\ast ,1)(i,1,m+1)}& {\alpha }_2{\delta }_{(i,\ast ,1)(i,1,,m+1)}& {\alpha }_3{\delta }_{(i,\ast ,1)(i,1,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(i,\ast ,1)(i,1,m+1)}\\ {\alpha }_1{\delta }_{(i,\ast ,2)(i,1,m+1)}& {\alpha }_2{\delta }_{(i,\ast ,2)(i,1,,m+1)}& {\alpha }_3{\delta }_{(i,\ast ,2)(i,1,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(i,\ast ,2)(i,1,m+1)}\\ {\alpha }_1{\delta }_{(i,\ast ,3)(i,1,m+1)}& {\alpha }_2{\delta }_{(i,\ast ,3)(i,1,,m+1)}& {\alpha }_3{\delta }_{(i,\ast ,3)(i,1,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(i,\ast ,3)(i,1,m+1)}\\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ {\alpha }_1{\delta }_{(i,\ast ,m)(i,1,m+1)}& {\alpha }_2{\delta }_{(i,\ast ,m)(i,1,m+1)}& {\alpha }_3{\delta }_{(i,\ast ,m)(i,1,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(i,\ast ,m)(i,1,m+1)}\end{array}\right]$$

where i varies from 1 to $n_1$. $K_i$ is a square matrix of order m.

Let

$$L_j=\left[\begin{array}{ccccc}{\alpha }_1{\delta }_{(\ast ,j,1)(1,j,m+1)}& {\alpha }_2{\delta }_{(\ast ,j,1)(1,j,m+1)}& {\alpha }_3{\delta }_{(\ast ,j,1)(1,j,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(\ast ,j,1)(1,j,m+1)}\\ {\alpha }_1{\delta }_{(\ast ,j,2)(i,j,m+1)}& {\alpha }_2{\delta }_{(\ast ,j,2)(1,j,,m+1)}& {\alpha }_3{\delta }_{(\ast ,j,2)(1,j,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(\ast ,j,2)(1,j,m+1)}\\ {\alpha }_1{\delta }_{(\ast ,j,3)(\ast ,j,m+1)}& {\alpha }_2{\delta }_{(\ast ,j,3)(\ast ,j,m+1)}& {\alpha }_3{\delta }_{(\ast ,j,3)(1,j,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(\ast ,j,3)(1,j,m+1)}\\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ {\alpha }_1{\delta }_{(\ast ,j,m)(1,j,m+1)}& {\alpha }_2{\delta }_{(\ast ,j,m)(1,j,m+1)}& {\alpha }_3{\delta }_{(\ast ,j,m)(1,j,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(\ast ,j,m)(1,j,m+1)}\end{array}\right]$$

where j varies from 1 to $n_2$. $L_j$ is a square matrix of order m.

${\mathcal{K}}_i=\left[\begin{array}{ccccc}{K}_i& \mathbf{0}& \mathbf{0}& \cdots \cdots & \mathbf{0}\end{array}\right]$ is a matrix of order $m\times n_{2}m$.

$\mathcal{L}=\left[\begin{array}{ccccc}{L}_1& \mathbf{0}& \mathbf{0}& \cdots \cdots & \mathbf{0}\\ \mathbf{0}& L_2& \mathbf{0}& \cdots \cdots & \mathbf{0}\\ \mathbf{0}& \mathbf{0}& L_3& \cdots \cdots & \mathbf{0}\\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \cdots \cdots & L_{n_2}\end{array}\right]$ is a square matrix of order $n_{2}m$.

With these notations, ${\mathcal{C}}_0$ takes the form

$${\mathcal{C}}_0=\left[\begin{array}{ccccc}{\mathcal{K}}_1& \mathbf{0}& \mathbf{0}& \cdots \cdots & \mathbf{0}\\ \mathbf{0}& {\mathcal{K}}_2& \mathbf{0}& \cdots \cdots & \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathcal{K}}_3& \cdots \cdots & \mathbf{0}\\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \cdots \cdots & {\mathcal{K}}_{n_1}\\ \mathcal{L}& \mathbf{0}& \mathbf{0}& \cdots \cdots & \mathbf{0}\end{array}\right].$$

Let

$$M_i=\left[\begin{array}{cccc}{\delta }_{(i,n_2,1)(i,n_2+1,1)}& & & \\ & {\delta }_{(i,n_2,2)(i,n_2+1,2)}& & & \\ & & {\delta }_{(i,n_2,3)(i,n_2+1,3)}& & \\ & & & \ddots \\ & & & & {\delta }_{(i,n_2,m)(i,n_2+1,m)}\end{array}\right]$$

be a square matrix of order m for $1\le i\le n_1$. Furthermore, let $P_i=\left[\begin{array}{c}\mathbf{0}\\ \mathbf{0}\\ ⋮\\ M_i\end{array}\right]$ be a matrix of order $n_{2}m\times m$ with $1\le i\le n_1$.

We define

$$Q_j=\left[\begin{array}{cccc}{\delta }_{(n_1,j,1)(n_1+1,j,1)}& & & \\ & {\delta }_{(n_1,j,2)(n_1+1,j,2)}& & & \\ & & {\delta }_{(n_1,j,3)(n_1+1,j,3)}& & \\ & & & \ddots \\ & & & & {\delta }_{(n_1,j,m)(n_1+1,j,m)}\end{array}\right]$$

as a square matrix of order m for $1\le j\le n_2$.

$V=\left[\begin{array}{ccccc}{Q}_1& \mathbf{0}& \mathbf{0}& \cdots \cdots & \mathbf{0}\\ \mathbf{0}& Q_2& \mathbf{0}& \cdots \cdots & \mathbf{0}\\ \mathbf{0}& \mathbf{0}& Q_3& \cdots \cdots & \mathbf{0}\\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \cdots \cdots & Q_{n_2}\end{array}\right]$ is a square matrix of order $n_{2}m$.

Then,

$$C_2=\left[\begin{array}{cccccccc}{P}_1& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \cdots \cdots & \mathbf{0}& \mathbf{0}\\ \mathbf{0}& P_2& \mathbf{0}& \mathbf{0}& \mathbf{0}& \cdots \cdots & \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& P_3& \mathbf{0}& \mathbf{0}& \cdots \cdots & \mathbf{0}& \mathbf{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 \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \cdots \cdots & P_{n_1}& V\end{array}\right].$$

We define

$$S_{ij}=\left[\begin{array}{ccccc}{S}_{ij1}& {\delta }_{(i,j,1)(i,j,2)}& {\delta }_{(i,j,1)(i,j,3)}& \cdots \cdots & {\delta }_{(i,j,1)(i,j,m)}\\ {\delta }_{(i,j,2)(i,j,1)}& S_{ij2}& {\delta }_{(i,j,2)(i,j,,3)}& \cdots \cdots & {\delta }_{(i,j,2)(i,j,,m)}\\ {\delta }_{(i,j,3)(i,j,1)}& {\delta }_{(i,j,3)(i,j,2)}& S_{ij3}& \cdots \cdots & {\delta }_{(i,j,3)(i,j,m)}\\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ {\delta }_{(i,j,m)(i,j,1)}& {\delta }_{(i,j,,m)(i,j,2)}& {\delta }_{(i,j,m)(i,j,3)}& \cdots \cdots & S_{ijm}\end{array}\right]$$

as a square matrix of order m with $1\le j\le n_2$.

In the above, $S_{ijh}$ = $-[{\displaystyle \sum _{k=1,h\ne k}^{m+1}}{\delta }_{(i,j,h)(i,j,k)}+{\delta }_{(i,j,h)(i,j+1,h)}+{\delta }_{(i,j,h)(i+1,j,h)}]$ for $1\le h\le m$ with $1\le j\le n_2$.

Take

$$T_{ij}=\left[\begin{array}{cccc}{\delta }_{(i,j,1)(i,j+1,1)}& & & \\ & {\delta }_{(i,j,2)(i,j+1,2)}& & & \\ & & {\delta }_{(i,j,3)(i,j+1,3)}& & \\ & & & \ddots \\ & & & & {\delta }_{(i,j,m)(i,j+1,m)}\end{array}\right],$$

which is a square matrix of order m with $1\le j\le n_2-1.$

Let $N_i=\left[\begin{array}{ccccccc}{S}_{i1}& T_{i1}& & & & & \\ & S_{i2}& T_{i2}& & & & \\ & & S_{i3}& T_{i3}& & & & \\ & & & \ddots \ddots & \\ & & & & & S_{i(n_2-1}& T_{i(n_2-1)}& \\ & & & & & & S_{i{n}_2}\end{array}\right]$ be a square matrix of order $m{n}_2$ with $1\le i\le n_1$.

Let $U_{ij}$ $=\left[\begin{array}{cccc}{\delta }_{(i,j,1)(i+1,j,1)}& & & \\ & {\delta }_{(i,j,2)(i+1,j,2)}& & & \\ & & {\delta }_{(i,j,3)(i+1,j,3)}& & \\ & & & \ddots \\ & & & & {\delta }_{(i,j,m)(i+1,j,m)}\end{array}\right]$ be a square matrix of order m for $1\le j\le n_2.$

Let $W_i=\left[\begin{array}{ccccccc}{U}_{i1}& & & & & & \\ & U_{i2}& & & & & \\ & & U_{i3}& & & & & \\ & & & \ddots \ddots & \\ & & & & & U_{i(n_2-1)}& & \\ & & & & & & U_{i{n}_2}\end{array}\right]$ be a square matrix of order $m{n}_2$ with $1\le i\le n_1-1$.

Using the above notations, ${\mathcal{A}}_1$ can be written as

$${\mathcal{A}}_1=\left[\begin{array}{ccccccc}{N}_1& W_1& & & & & \\ & N_2& W_2& & & & \\ & & N_3& W_3& & & & \\ & & & \ddots \ddots & \\ & & & & & N_{n_1-1}& W_{n_1-1}& \\ & & & & & & N_{n_1}\end{array}\right].$$

Let ${\mathcal{D}}_i=\left[\begin{array}{ccccc}{P}_i& \mathbf{0}& \mathbf{0}& \cdots \cdots & \mathbf{0}\end{array}\right]$ be a square matrix of order $n_{2}m$ for $1\le i\le n_1$.

Using these, ${\mathcal{A}}_2$ can be expressed as

$${\mathcal{A}}_2=\left[\begin{array}{ccccccc}{\mathcal{D}}_1& & & & & & \\ & {\mathcal{D}}_2& & & & & \\ & & {\mathcal{D}}_3& & & & & \\ & & & \ddots & \\ & & & & & {\mathcal{D}}_{n_1-1}& & \\ V& & & & & & {\mathcal{D}}_{n_1}\end{array}\right].$$

We define ${\mathcal{F}}_{ij}$ as follows:

$${\mathcal{F}}_{ij}=\left[\begin{array}{ccccc}{\alpha }_1{\delta }_{(i,j,1)(i,j,m+1)}& {\alpha }_2{\delta }_{(i,j,1)(i,j,m+1)}& {\alpha }_3{\delta }_{(i,j,1)(i,j,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(i,j,1)(i,j,m+1)}\\ {\alpha }_1{\delta }_{(i,j,2)(i,j,m+1)}& {\alpha }_2{\delta }_{(i,j,2)(i,j,m+1)}& {\alpha }_3{\delta }_{(i,j,2)(i,j,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(i,j,2)(i,j,m+1)}\\ {\alpha }_1{\delta }_{(i,j,3)(i,j,m+1)}& {\alpha }_2{\delta }_{(i,j,3)(i,j,m+1)}& {\alpha }_3{\delta }_{(i,j,3)(i,j,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(i,j,3)(i,j,m+1)}\\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots & \cdots \cdots \\ {\alpha }_1{\delta }_{(i,j,m)(i,j,m+1)}& {\alpha }_2{\delta }_{(i,j,m)(i,j,m+1)}& {\alpha }_3{\delta }_{(i,j,m)(i,j,m+1)}& \cdots \cdots & {\alpha }_m{\delta }_{(i,j,m)(i,j,m+1)}\end{array}\right]$$

where j varies from 1 to $n_2$. ${\mathcal{F}}_{ij}$ is a square matrix of order m.

Then, we can write

$${\mathcal{G}}_i=\left[\begin{array}{ccccccc}{\mathcal{F}}_{i1}& & & & & & \\ & {\mathcal{F}}_{i2}& & & & & \\ & & {\mathcal{F}}_{i3}& & & & & \\ & & & \ddots \ddots & \\ & & & & & {\mathcal{F}}_{i{n}_2-1}& & \\ & & & & & & {\mathcal{F}}_{i{n}_2}\end{array}\right],$$

which is a square matrix of order $m{n}_2$ with $1\le i\le n_1$.

These help us to represent ${\mathcal{A}}_0$ as follows:

$${\mathcal{A}}_0=\left[\begin{array}{ccccccc}{\mathcal{G}}_1& & & & & & \\ & {\mathcal{G}}_2& & & & & \\ & & {\mathcal{G}}_3& & & & & \\ & & & \ddots & \\ & & & & & {\mathcal{G}}_{n_1-1}& & \\ & & & & & & {\mathcal{G}}_{n_1}\end{array}\right].$$

Next, we proceed to compute the stability condition of the system and obtain the system state probability vectors. From each of these vectors, the system occupying a specific state can be computed.

3. Steady State Analysis

In this section, we perform steady-state analysis of the queueing model by initially determining the stability condition of the system.

3.1. Stability Condition

The infinitesimal generator $\mathcal{A}={\mathcal{A}}_0+{\mathcal{A}}_1+{\mathcal{A}}_2$ is

$$\left[\begin{array}{ccccccc}{\mathcal{G}}_1+{\mathcal{N}}_1+{\mathcal{D}}_1& {\mathcal{U}}_1& & & & & \\ & {\mathcal{G}}_2+{\mathcal{N}}_2+{\mathcal{D}}_2& {\mathcal{U}}_2& & & & \\ & & {\mathcal{G}}_3+{\mathcal{N}}_3+{\mathcal{D}}_3& {\mathcal{U}}_3& & & \\ & & & \ddots \ddots & \\ & & & & & {\mathcal{G}}_{n_1-1}+{\mathcal{N}}_{n_1-1}+{\mathcal{D}}_{n_1-1}& {\mathcal{U}}_{n_1}-1\\ V& & & & & & {\mathcal{G}}_{n_1}+{\mathcal{N}}_{n_1}+{\mathcal{D}}_{n_1},\end{array}\right]$$

which is a square matrix of order $n_1{n}_{2}m$

Let $\mathit{\pi }=({\mathit{\pi }}_1,{\mathit{\pi }}_2,{\mathit{\pi }}_3,\dots \dots ,{\mathit{\pi }}_{n_1})$ denote the steady-state probability vector of the generator matrix $\mathcal{A}$.

Here, $\mathit{\pi }$ is of dimension $1\times n_1{n}_{2}m$ and ${\mathit{\pi }}_r$ is of dimension $1\times n_{2}m$ for $r=1,2,\dots n_1$.

The steady state probability vector $\mathit{\pi }$ satisfies the equations

$$\mathit{\pi }\mathcal{A}=0,\mathit{\pi }\mathbf{e}=1.$$

(1)

Using Equation (1), we obtain

$$\begin{array}{cc}\hfill {\mathit{\pi }}_1[{\mathcal{G}}_1+{\mathcal{N}}_1+{\mathcal{D}}_1]+{\mathit{\pi }}_{n_1}V& \hfill =\mathbf{0}\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {\mathit{\pi }}_1{\mathcal{U}}_1+{\mathit{\pi }}_2[{\mathcal{G}}_2+{\mathcal{N}}_2+{\mathcal{D}}_2]& \hfill =\mathbf{0}\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\mathit{\pi }}_2{\mathcal{U}}_2+{\mathit{\pi }}_3[{\mathcal{G}}_3+{\mathcal{N}}_3+{\mathcal{D}}_3]& \hfill =\mathbf{0}\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\mathit{\pi }}_{n_1-1}{\mathcal{U}}_{n_1-1}+{\mathit{\pi }}_{n_1}[{\mathcal{G}}_{n_1}+{\mathcal{N}}_{n_1}+{\mathcal{D}}_{n_1}]& \hfill =\mathbf{0}\hfill \end{array}$$

(5)

$$\begin{array}{c}\hfill {\mathit{\pi }}_1\times \mathbf{e}+{\mathit{\pi }}_2\times \mathbf{e}+\dots .+{\mathit{\pi }}_{n_1}\times \mathbf{e}=1.\end{array}$$

(6)

By solving the above system equations, we can find ${\mathit{\pi }}_1,{\mathit{\pi }}_2,{\mathit{\pi }}_3,\dots \dots {\mathit{\pi }}_{n_1}$.

The LIQBD description of the model indicates that the queueing system is stable if

and only if the left drift exceeds the right drift (see Theorem 3.1.1 of Neuts [35]), that is,

$$\mathit{\pi }{\mathcal{A}}_0\mathbf{e}<\mathit{\pi }{\mathcal{A}}_2\mathbf{e},$$

(7)

$$\mathit{\pi }{\mathcal{A}}_0\mathbf{e}={\mathit{\pi }}_1{\mathcal{G}}_1\mathbf{e}+{\mathit{\pi }}_2{\mathcal{G}}_2\mathbf{e}+{\mathit{\pi }}_3{\mathcal{G}}_3\mathbf{e}+\dots +{\mathit{\pi }}_{n_1}{\mathcal{G}}_{n_1}\mathbf{e}=\sum _{i=1}^{n_1}{\mathit{\pi }}_i{\mathcal{G}}_i\mathbf{e},$$

(8)

$$\begin{array}{c}\hfill \mathit{\pi }{A}_2\mathbf{e}={\mathit{\pi }}_1{\mathcal{D}}_1\mathbf{e}+{\mathit{\pi }}_2{\mathcal{D}}_2\mathbf{e}+{\mathit{\pi }}_3{\mathcal{D}}_3\mathbf{e}+\dots +{\mathit{\pi }}_{n_1}{\mathcal{D}}_{n_1}\mathbf{e}+{\mathit{\pi }}_{n_1}V\mathbf{e}=\sum _{i=1}^{n_1}{\mathit{\pi }}_i{\mathcal{D}}_i\mathbf{e}+{\mathit{\pi }}_{n_1}V\mathbf{e}.\end{array}$$

(9)

Theorem 1. 

The given system is stable if and only if

$$\sum _{i=1}^{n_1}{\mathit{\pi }}_i{\mathcal{G}}_i\mathbf{e}<\sum _{i=1}^{n_1}{\mathit{\pi }}_i{\mathcal{D}}_i\mathbf{e}+{\mathit{\pi }}_{n_1}V\mathbf{e}.$$

(10)

3.2. The Steady State Probability Vector of $Q^{\ast }$

Let $\mathit{x}$ be the steady state probability vector of $Q^{\ast }$: $\mathit{x}=({\mathit{x}}_0,{\mathit{x}}_1,{\mathit{x}}_2\dots )$, where ${\mathit{x}}_0$ is of dimension $1\times m$, ${\mathit{x}}_1$ is of dimension $1\times (n_1+n_2)m$, and ${\mathit{x}}_2,{\mathit{x}}_3,\dots {\mathit{x}}_{n_1}$ are of dimension $1\times n_1{n}_{2}m$.

Under the stability condition, we have ${\mathit{x}}_{\mathit{i}}={\mathit{x}}_2{R}^{i-2},i\ge 3$, where the matrix R is the minimal nonnegative solution to the matrix quadratic equation

$$R^2{\mathcal{A}}_2+R{\mathcal{A}}_1+{\mathcal{A}}_0=0$$

and the vectors ${\mathit{x}}_0$, ${\mathit{x}}_1$ and ${\mathit{x}}_2$ are obtained by solving the equations

$$\begin{array}{c}\hfill {\mathit{x}}_0{\mathcal{B}}_1+{\mathit{x}}_1{\mathcal{B}}_2=0\end{array}$$

(11)

$$\begin{array}{c}\hfill {\mathit{x}}_0{\mathcal{B}}_0+{\mathit{x}}_1{\mathcal{C}}_1+{\mathit{x}}_{\mathbf{2}}{\mathcal{C}}_2=0\end{array}$$

(12)

$$\begin{array}{c}\hfill {\mathit{x}}_1{\mathcal{C}}_0+{\mathit{x}}_2({\mathcal{A}}_1+R{\mathcal{A}}_2)=0\end{array}$$

(13)

subject to the normalizing condition

$${\mathit{x}}_0\mathit{e}+{\mathit{x}}_1\mathit{e}+{\mathit{x}}_2{(I-R)}^{-1}\mathit{e}=1.$$

(14)

Having obtained the system state probability vectors, we can now make use of them to compute important performance measures.

4. Performance Measures

Formulae for key performance measures are derived.

• The probability that both of the servers are idle: $p_{idle}={\mathit{x}}_0\mathbf{e}$.
• The probability that only Server 1 ($S_1$) is idle:

  $$p_{1idle}=\sum _{j=1}^{n_2}\sum _{k=1}^m{\mathit{x}}_{1\ast jk}$$
• The probability that only Server 2 ($S_2$) is idle:

  $$p_{2idle}=\sum _{i=1}^{n_1}\sum _{k=1}^m{\mathit{x}}_{1i\ast k}$$
• The probability that Server 1 is busy:

  $$P_{1busy}=\sum _{q=2}^{\infty }\sum _{i=1}^{n_1}\sum _{j=1}^{n_2}\sum _{k=1}^m{\mathit{x}}_{qijk}+\sum _{i=1}^{n_1}\sum _{k=1}^m{\mathit{x}}_{1i\ast k}$$
• The probability that Server 2 is busy:

  $$P_{2busy}=\sum _{q=2}^{\infty }\sum _{i=1}^{n_1}\sum _{j=1}^{n_2}\sum _{k=1}^m{\mathit{x}}_{qijk}+\sum _{j=1}^{n_2}\sum _{k=1}^m{\mathit{x}}_{1\ast jk}$$
• The probability that both the servers are busy:

  $$P_{busy}=\sum _{q=2}^{\infty }\sum _{i=1}^{n_1}\sum _{j=1}^{n_2}\sum _{k=1}^m{\mathit{x}}_{qijk}$$
• The probability of q customers being in the system: $P_q={\mathit{x}}_q\mathbf{e}$.
• The expected number of customers in the system:

  $$ECS=\sum _{q=1}^{\infty }q{\mathit{x}}_q\mathbf{e}$$
• The expected number of customers in the queue:
  ${\displaystyle ECQ=\sum _{q=2}^{\infty }(q-2){\mathit{x}}_q\mathbf{e}}$

5. Numerical Illustration

In this section, we provide four numerical illustrations of our model.

We take $n_1=2,n_2=3,$ and $m=2$. The state space of the arrival process is $1,2,3$, where the transient states are 1,2 and the absorbing state is 3. The state space of the service process provided by server 1 ($S_1$) is $\{1,2,3\}$, where the transient states are 1,2 and the absorbing state is 3. The state space of the service process provided by server 2 ($S_2$) is $\{1,2,3,4\}$, where the transient states are 1,2,3 and the absorbing state is 4. Then, the state space of the Markov chain $\mathcal{X}$ = $\left\{\right(i,j,k):1\le i\le 3;1\le j\le 4;1\le k\le 3\}$. The absorbing states of $\mathcal{X}$ are: $\left\{\right\{(i,j,3):1\le i\le 3;1\le j\le 4\}\bigcup \{(i,4,k):1\le i\le 3;1\le k\le 3\}\bigcup \{(3,j,k):1\le j\le 4;$ $1\le k\le 3\left\}\right\}\bigcup \left\{\right(3,4,k):1\le k\le 3\}\bigcup \left\{\right(3,j,3):1\le j\le 4\}\bigcup \left\{\right(i,4,3):1\le i\le 3\}\bigcup \left\{\right(3,4,3\left)\right\}$.

We assume that at most one coordinate change in a transition has a positive probability. In the absence of customers, no service can be provided, as indicated by a * symbol in the position of the service coordinates. Here, $\mathit{\alpha }$ = $(0.7,0.3)$, ${\alpha }_{m+1}=0$. The case of more than one coordinate change in a single transition (only up to three being possible in the present case, as we have the service stages provided by each server as first two coordinates and the changes in the arrival phases as the last coordinate) can be considered in a similar fashion. Consequently, we obtain a much less sparse transition rate matrix.

5.1. Illustration 1

The transition rates are taken as shown in

Table 2. Transition rates.

	$(2,1,1)$	$(2,1,2)$	$(2,2,1)$	$(2,2,2)$	$(2,3,1)$	$(2,3,2)$	$(2,1,3)$	$(2,2,3)$	$(2,3,3)$	$(2,4,1)$	$(2,4,2)$	$(3,1,1)$	$(3,1,2)$	$(3,2,1)$	$(3,2,2)$	$(3,3,1)$	$(3,3,2)$
$(1,1,1)$	3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,1,2)$	0	3.1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,2,1)$	0	0	3.1	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,2,2)$	0	0	0	2.3	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,3,1)$	0	0	0	0	2.8	0	0	0	0	0	0	0	0	0	0	0	0
$(1,3,2)$	0	0	0	0	0	2.3	0	0	0	0	0	0	0	0	0	0	0
$(2,1,1)$	−10.7	3.1	2.8	0	0	0	1.7	0	0	0	0	3.1	0	0	0	0	0
$(2,1,2)$	2.6	−12.3	0	3.1	0	0	3.1	0	0	0	0	0	3.5	0	0	0	0
$(2,2,1)$	0	0	−13.1	2.4	3.8	0	0	3.2	0	0	0	0	0	3.7	0	0	0
$(2,2,2)$	0	0	2.7	−11.1	0	2.2	0	3.1	0	0	0	0	0	0	3.1	0	0
$(2,3,1)$	0	0	0	0	−9.5	1.1	0	0	2.6	3	0	0	0	0	0	2.8	0
$(2,3,2)$	0	0	0	0	1.9	−12.1	0	0	2.5	0	4.2	0	0	0	0	0	3.5

and

Table 3. Transitions rates.

	$(\ast ,\ast ,1)$	$(\ast ,\ast ,2)$	$(1,\ast ,1)$	$(1,\ast ,2)$	$(2,\ast ,1)$	$(2,\ast ,2)$	$(\ast ,1,1)$	$(\ast ,1,2)$	$(\ast ,2,1)$	$(\ast ,2,2)$	$(\ast ,3,1)$	$(\ast ,3,2)$	$(1,1,1)$	$(1,1,2)$	$(1,2,1)$	$(1,2,2)$	$(1,3,1)$	$(1,3,2)$	$(1,4,1)$	$(1,4,2)$	$(1,\ast ,3)$	$(3,\ast ,1)$	$(3,\ast ,2)$	$(\ast ,4,1)$	$(\ast ,4,2)$	$(1,1,3)$	$(1,2,3)$	$(1,3,3)$	$(2,1,3)$
$(\ast ,\ast ,1)$	−4.5	1.3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	3.2	0	0	0	0	0	0	0	0
$(\ast ,\ast ,2)$	2.4	−5.2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.8	0	0	0	0	0	0	0	0
$(1,\ast ,1)$	0	0	−8.1	3.2	3.1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.8	0	0	0
$(1,\ast ,2)$	0	0	2.4	−6.7	0	2.7	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.6	0	0	0
$(2,\ast ,1)$	0	0	0	0	−6.5	1.5	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.8	0	0	0	0	0	0	2.2
$(2,\ast ,2)$	0	0	0	0	2.1	−5.9	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.9	0	0	0	0	0	1.9
$(\ast ,1,1)$	0	0	0	0	0	0	−6.2	1.7	2.6	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.9	0	0	0
$(\ast ,1,2)$	0	0	0	0	0	0	2.2	−6.4	0	2.4	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.8	0	0	0
$(\ast ,2,1)$	0	0	0	0	0	0	0	0	−6.1	2.4	2.1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.6	0	0
$(\ast ,2,2)$	0	0	0	0	0	0	0	0	2.9	−8	0	2.9	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.2	0	0
$(\ast ,3,1)$	0	0	0	0	0	0	0	0	0	0	−6.7	2.2	0	0	0	0	0	0	0	0	0	0	0	2.1	0	0	0	2.4	0
$(\ast ,3,2)$	0	0	0	0	0	0	0	0	0	0	2.1	−8.5	0	0	0	0	0	0	0	0	0	0	0	0	3.8	0	0	2.6	0
$(1,1,1)$	0	0	0	0	0	0	0	0	0	0	0	0	−11.2	3	2.6	0	0	0	0	0	0	0	0	0	0	2.6	0	0	0
$(1,1,2)$	0	0	0	0	0	0	0	0	0	0	0	0	2.8	−11.7	0	2.7	0	0	0	0	0	0	0	0	0	3.1	0	0	0
$(1,2,1)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	−8.9	1.5	2.2	0	0	0	0	0	0	0	0	0	2.1	0	0
$(1,2,2)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	3.2	−10	0	2.5	0	0	0	0	0	0	0	0	2	0	0
$(1,3,1)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	−10.7	1.9	4.2	0	0	0	0	0	0	0	0	1.8	0
$(1,3,2)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.5	−11.4	0	3.4	0	0	0	0	0	0	0	3.2	0

.

Using the transition rates shown in Table 2 and Table 3, we obtain the values of the performance measures as follows.

$$\begin{array}{cc}\hfill P_{idle}& =0.0134\hfill \\ \hfill P_{1idle}& =0.0420\hfill \\ \hfill P_{2idle}& =0.0365\hfill \\ \hfill P_{1busy}& =0.9445\hfill \\ \hfill P_{2busy}& =0.9501\hfill \\ \hfill P_{busy}& =0.9081\hfill \\ \hfill ECS& =10.7079\hfill \\ \hfill ECQ& =8.8133\hfill \end{array}$$

5.2. Illustration 2

The transition rates are taken as shown in

Table 4. Transition rates.

	$(2,1,1)$	$(2,1,2)$	$(2,2,1)$	$(2,2,2)$	$(2,3,1)$	$(2,3,2)$	$(2,1,3)$	$(2,2,3)$	$(2,3,3)$	$(2,4,1)$	$(2,4,2)$	$(3,1,1)$	$(3,1,2)$	$(3,2,1)$	$(3,2,2)$	$(3,3,1)$	$(3,3,2)$
$(1,1,1)$	3.2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,1,2)$	0	3.2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,2,1)$	0	0	3.3	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,2,2)$	0	0	0	2.5	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,3,1)$	0	0	0	0	2.7	0	0	0	0	0	0	0	0	0	0	0	0
$(1,3,2)$	0	0	0	0	0	2.6	0	0	0	0	0	0	0	0	0	0	0
$(2,1,1)$	−11.1	3.5	2.6	0	0	0	1.8	0	0	0	0	3.2	0	0	0	0	0
$(2,1,2)$	2.4	−11.9	0	3.3	0	0	3	0	0	0	0	0	3.2	0	0	0	0
$(2,2,1)$	0	0	−13.2	2.6	3.2	0	0	3.5	0	0	0	0	0	3.9	0	0	0
$(2,2,2)$	0	0	2.5	−12.8	0	3.4	0	3.5	0	0	0	0	0	0	3.4	0	0
$(2,3,1)$	0	0	0	0	−9.3	1	0	0	2.9	3.1	0	0	0	0	0	2.3	0
$(2,3,2)$	0	0	0	0	1.5	−12.3	0	0	2.2	0	4.9	0	0	0	0	0	3.7

and

Table 5. Transition rates.

	$(\ast ,\ast ,1)$	$(\ast ,\ast ,2)$	$(1,\ast ,1)$	$(1,\ast ,2)$	$(2,\ast ,1)$	$(2,\ast ,2)$	$(\ast ,1,1)$	$(\ast ,1,2)$	$(\ast ,2,1)$	$(\ast ,2,2)$	$(\ast ,3,1)$	$(\ast ,3,2)$	$(1,1,1)$	$(1,1,2)$	$(1,2,1)$	$(1,2,2)$	$(1,3,1)$	$(1,3,2)$	$(1,4,1)$	$(1,4,2)$	$(1,\ast ,3)$	$(3,\ast ,1)$	$(3,\ast ,2)$	$(\ast ,4,1)$	$(\ast ,4,2)$	$(1,1,3)$	$(1,2,3)$	$(1,3,3)$	$(2,1,3)$
$(\ast ,\ast ,1)$	−4.6	1.5	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	3.1	0	0	0	0	0	0	0	0
$(\ast ,\ast ,2)$	2	−4.5	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.5	0	0	0	0	0	0	0	0
$(1,\ast ,1)$	0	0	−8.9	3.6	3.4	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.9	0	0	0
$(1,\ast ,2)$	0	0	2.9	−7.3	0	2.9	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.5	0	0	0
$(2,\ast ,1)$	0	0	0	0	−6.1	2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0	0	0	0	2.1
$(2,\ast ,2)$	0	0	0	0	2.8	−5.9	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.5	0	0	0	0	0	1.6
$(\ast ,1,1)$	0	0	0	0	0	0	−6.1	1.9	2.1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.1	0	0	0
$(\ast ,1,2)$	0	0	0	0	0	0	2.7	−6.8	0	2.6	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.5	0	0	0
$(\ast ,2,1)$	0	0	0	0	0	0	0	0	−6.2	2.5	2.3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.4	0	0
$(\ast ,2,2)$	0	0	0	0	0	0	0	0	2.7	−7.4	0	2.6	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.1	0	0
$(\ast ,3,1)$	0	0	0	0	0	0	0	0	0	0	−7	2.4	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0	2.6	0
$(\ast ,3,2)$	0	0	0	0	0	0	0	0	0	0	1.9	−8.6	0	0	0	0	0	0	0	0	0	0	0	0	3.9	0	0	2.8	0
$(1,1,1)$	0	0	0	0	0	0	0	0	0	0	0	0	−10.2	3.1	2.3	0	0	0	0	0	0	0	0	0	0	1.6	0	0	0
$(1,1,2)$	0	0	0	0	0	0	0	0	0	0	0	0	2.5	−11.1	0	2.9	0	0	0	0	0	0	0	0	0	2.5	0	0	0
$(1,2,1)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	−10.5	1.7	3.2	0	0	0	0	0	0	0	0	0	2.3	0	0
$(1,2,2)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.2	−8.5	0	2.1	0	0	0	0	0	0	0	0	1.7	0	0
$(1,3,1)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	−12.1	2.1	4.5	0	0	0	0	0	0	0	0	2.8	0
$(1,3,2)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.8	−11.6	0	3.2	0	0	0	0	0	0	0	3	0

.

Using the transition rates shown in Table 4 and Table 5, we obtain the values of the performance measures are as follows.

$$\begin{array}{cc}\hfill P_{idle}& =0.0040\hfill \\ \hfill P_{1idle}& =0.0122\hfill \\ \hfill P_{2idle}& =0.0113\hfill \\ \hfill P_{1busy}& =0.9838\hfill \\ \hfill P_{2busy}& =0.9847\hfill \\ \hfill P_{busy}& =0.9725\hfill \\ \hfill ECS& =40.4592\hfill \\ \hfill ECQ& =38.4907\hfill \end{array}$$

5.3. Illustration 3

The transition Rates are taken as shown in

Table 6. Transition rates.

	$(2,1,1)$	$(2,1,2)$	$(2,2,1)$	$(2,2,2)$	$(2,3,1)$	$(2,3,2)$	$(2,1,3)$	$(2,2,3)$	$(2,3,3)$	$(2,4,1)$	$(2,4,2)$	$(3,1,1)$	$(3,1,2)$	$(3,2,1)$	$(3,2,2)$	$(3,3,1)$	$(3,3,2)$
$(1,1,1)$	3.2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,1,2)$	0	3.7	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,2,1)$	0	0	3.5	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,2,2)$	0	0	0	2.7	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,3,1)$	0	0	0	0	2.9	0	0	0	0	0	0	0	0	0	0	0	0
$(1,3,2)$	0	0	0	0	0	2.9	0	0	0	0	0	0	0	0	0	0	0
$(2,1,1)$	−12.2	3.8	2.9	0	0	0	2.1	0	0	0	0	3.4	0	0	0	0	0
$(2,1,2)$	2.7	−12.8	0	3.3	0	0	3.2	0	0	0	0	0	3.6	0	0	0	0
$(2,2,1)$	0	0	−15.4	2.9	3.8	0	0	3.8	0	0	0	0	0	4.9	0	0	0
$(2,2,2)$	0	0	2.8	−12.6	0	3	0	3.6	0	0	0	0	0	0	3.2	0	0
$(2,3,1)$	0	0	0	0	−11.2	1.3	0	0	3.1	3.9	0	0	0	0	0	2.9	0
$(2,3,2)$	0	0	0	0	2	−13.9	0	0	2.6	0	5.4	0	0	0	0	0	3.9

and

Table 7. Transitions rates.

	$(\ast ,\ast ,1)$	$(\ast ,\ast ,2)$	$(1,\ast ,1)$	$(1,\ast ,2)$	$(2,\ast ,1)$	$(2,\ast ,2)$	$(\ast ,1,1)$	$(\ast ,1,2)$	$(\ast ,2,1)$	$(\ast ,2,2)$	$(\ast ,3,1)$	$(\ast ,3,2)$	$(1,1,1)$	$(1,1,2)$	$(1,2,1)$	$(1,2,2)$	$(1,3,1)$	$(1,3,2)$	$(1,4,1)$	$(1,4,2)$	$(1,\ast ,3)$	$(3,\ast ,1)$	$(3,\ast ,2)$	$(\ast ,4,1)$	$(\ast ,4,2)$	$(1,1,3)$	$(1,2,3)$	$(1,3,3)$	$(2,1,3)$
$(\ast ,\ast ,1)$	−5.1	1.8	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	3.3	0	0	0	0	0	0	0	0
$(\ast ,\ast ,2)$	2.4	−5	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.6	0	0	0	0	0	0	0	0
$(1,\ast ,1)$	0	0	−9.2	3.8	3.5	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.9	0	0	0
$(1,\ast ,2)$	0	0	3.1	−7.9	0	3.2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.6	0	0	0
$(2,\ast ,1)$	0	0	0	0	−8	2.3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	3.5	0	0	0	0	0	0	2.2
$(2,\ast ,2)$	0	0	0	0	3	−6.6	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.9	0	0	0	0	0	1.7
$(\ast ,1,1)$	0	0	0	0	0	0	−6.8	2.2	2.4	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.2	0	0	0
$(\ast ,1,2)$	0	0	0	0	0	0	2.9	−7.2	0	2.7	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.6	0	0	0
$(\ast ,2,1)$	0	0	0	0	0	0	0	0	−6.6	2.6	2.4	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.6	0	0
$(\ast ,2,2)$	0	0	0	0	0	0	0	0	2.9	−7.9	0	2.8	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.2	0	0
$(\ast ,3,1)$	0	0	0	0	0	0	0	0	0	0	−7.7	2.6	0	0	0	0	0	0	0	0	0	0	0	2.2	0	0	0	2.9	0
$(\ast ,3,2)$	0	0	0	0	0	0	0	0	0	0	2.2	−9.4	0	0	0	0	0	0	0	0	0	0	0	0	4.1	0	0	3.1	0
$(1,1,1)$	0	0	0	0	0	0	0	0	0	0	0	0	−10.7	3.4	2.4	0	0	0	0	0	0	0	0	0	0	1.7	0	0	0
$(1,1,2)$	0	0	0	0	0	0	0	0	0	0	0	0	2.7	−12.3	0	3.1	0	0	0	0	0	0	0	0	0	2.8	0	0	0
$(1,2,1)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	−11.6	1.9	3.4	0	0	0	0	0	0	0	0	0	2.8	0	0
$(1,2,2)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.5	−9.5	0	2.4	0	0	0	0	0	0	0	0	1.9	0	0
$(1,3,1)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	−13.1	2.4	4.7	0	0	0	0	0	0	0	0	3.1	0
$(1,3,2)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	3	−12.7	0	3.5	0	0	0	0	0	0	0	3.3	0

.

Using the transition rates shown in Table 6 and Table 7, we obtain the values of the performance measures as follows.

$$\begin{array}{cc}\hfill P_{idle}& =0.0023\hfill \\ \hfill P_{1idle}& =0.0064\hfill \\ \hfill P_{2idle}& =0.0059\hfill \\ \hfill P_{1busy}& =0.9913\hfill \\ \hfill P_{2busy}& =0.9918\hfill \\ \hfill P_{busy}& =0.9854\hfill \\ \hfill ECS& =78.0583\hfill \\ \hfill ECQ& =76.0752\hfill \end{array}$$

5.4. Illustration 4

The transition rates are taken as shown in

Table 8. Transition rates.

	$(2,1,1)$	$(2,1,2)$	$(2,2,1)$	$(2,2,2)$	$(2,3,1)$	$(2,3,2)$	$(2,1,3)$	$(2,2,3)$	$(2,3,3)$	$(2,4,1)$	$(2,4,2)$	$(3,1,1)$	$(3,1,2)$	$(3,2,1)$	$(3,2,2)$	$(3,3,1)$	$(3,3,2)$
$(1,1,1)$	3.3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,1,2)$	0	3.8	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,2,1)$	0	0	3.6	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,2,2)$	0	0	0	2.8	0	0	0	0	0	0	0	0	0	0	0	0	0
$(1,3,1)$	0	0	0	0	3	0	0	0	0	0	0	0	0	0	0	0	0
$(1,3,2)$	0	0	0	0	0	3	0	0	0	0	0	0	0	0	0	0	0
$(2,1,1)$	−12.6	3.9	3	0	0	0	2.2	0	0	0	0	3.5	0	0	0	0	0
$(2,1,2)$	2.8	−13.2	0	3.4	0	0	3.3	0	0	0	0	0	3.7	0	0	0	0
$(2,2,1)$	0	0	−15.8	3	3.9	0	0	3.9	0	0	0	0	0	5	0	0	0
$(2,2,2)$	0	0	3	−13.4	0	3.2	0	3.8	0	0	0	0	0	0	3.4	0	0
$(2,3,1)$	0	0	0	0	−12.1	1.4	0	0	3.2	4.5	0	0	0	0	0	3	0
$(2,3,2)$	0	0	0	0	2.1	−14.7	0	0	2.7	0	5.9	0	0	0	0	0	4

and

Table 9. Transitions rates.

	$(\ast ,\ast ,1)$	$(\ast ,\ast ,2)$	$(1,\ast ,1)$	$(1,\ast ,2)$	$(2,\ast ,1)$	$(2,\ast ,2)$	$(\ast ,1,1)$	$(\ast ,1,2)$	$(\ast ,2,1)$	$(\ast ,2,2)$	$(\ast ,3,1)$	$(\ast ,3,2)$	$(1,1,1)$	$(1,1,2)$	$(1,2,1)$	$(1,2,2)$	$(1,3,1)$	$(1,3,2)$	$(1,4,1)$	$(1,4,2)$	$(1,\ast ,3)$	$(3,\ast ,1)$	$(3,\ast ,2)$	$(\ast ,4,1)$	$(\ast ,4,2)$	$(1,1,3)$	$(1,2,3)$	$(1,3,3)$	$(2,1,3)$
$(\ast ,\ast ,1)$	−5.3	1.9	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	3.4	0	0	0	0	0	0	0	0
$(\ast ,\ast ,2)$	2.5	−5.2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.7	0	0	0	0	0	0	0	0
$(1,\ast ,1)$	0	0	−9.5	3.9	3.6	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0
$(1,\ast ,2)$	0	0	3.2	−8.2	0	3.3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.7	0	0	0
$(2,\ast ,1)$	0	0	0	0	−8.3	2.4	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	3.6	0	0	0	0	0	0	2.3
$(2,\ast ,2)$	0	0	0	0	3.1	−6.9	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0	0	0	1.8
$(\ast ,1,1)$	0	0	0	0	0	0	−7.1	2.3	2.5	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.3	0	0	0
$(\ast ,1,2)$	0	0	0	0	0	0	3	−7.5	0	2.8	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.7	0	0	0
$(\ast ,2,1)$	0	0	0	0	0	0	0	0	−6.9	2.7	2.5	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1.7	0	0
$(\ast ,2,2)$	0	0	0	0	0	0	0	0	3	−8.2	0	2.9	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.3	0	0
$(\ast ,3,1)$	0	0	0	0	0	0	0	0	0	0	−8	2.7	0	0	0	0	0	0	0	0	0	0	0	2.3	0	0	0	3	0
$(\ast ,3,2)$	0	0	0	0	0	0	0	0	0	0	2.3	−9.7	0	0	0	0	0	0	0	0	0	0	0	0	4.2	0	0	3.2	0
$(1,1,1)$	0	0	0	0	0	0	0	0	0	0	0	0	−11.1	3.5	2.5	0	0	0	0	0	0	0	0	0	0	1.8	0	0	0
$(1,1,2)$	0	0	0	0	0	0	0	0	0	0	0	0	2.8	−12.7	0	3.2	0	0	0	0	0	0	0	0	0	2.9	0	0	0
$(1,2,1)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	−12	2	3.5	0	0	0	0	0	0	0	0	0	2.9	0	0
$(1,2,2)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2.6	−9.9	0	2.5	0	0	0	0	0	0	0	0	2	0	0
$(1,3,1)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	−13.5	2.5	4.8	0	0	0	0	0	0	0	0	3.2	0
$(1,3,2)$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	3.1	−13.1	0	3.6	0	0	0	0	0	0	0	3.4	0

.

Using the transition rates shown in Table 8 and Table 9, we obtain the values of the performance measures as follows.

$$\begin{array}{cc}\hfill P_{idle}& =0.0020\hfill \\ \hfill P_{1idle}& =0.0055\hfill \\ \hfill P_{2idle}& =0.0051\hfill \\ \hfill P_{1busy}& =0.9925\hfill \\ \hfill P_{2busy}& =0.9929\hfill \\ \hfill P_{busy}& =0.9874\hfill \\ \hfill ECS& =89.0997\hfill \\ \hfill ECQ& =87.1143\hfill \end{array}$$

From the above numerical illustrations, it can be concluded that when the transition rates increase, the values of $ECS$, $ECQ$, $P_{1busy}$, $P_{2busy}$, and $P_{busy}$ increase, while the values of $P_{idle}$, $P_{1idle}$, and $P_{2idle}$ decrease as the transition rates increases.

6. Conclusions

In this paper, we have analysed a queuing system in which the arrival process and services handled by two servers are interdependent. The evolution of these processes is governed by transitions occurring in the product space of three Markov chains according to a semi-Markov rule, with sojourn times in states following the exponential distribution. Using the matrix geometric method, we have analysed the model and computed various performance measures. Additionally, we have included numerical illustrations of the model.

We propose extending the model presented here to a more general case. In this paper, we have assumed that upon the arrival of a customer in the system, the starting phase of the next arrival is chosen according to the initial probability vector. It is possible to extend this to the case in which the next arrival phase is the phase in which transition occurs consequent to the arrival of the customer.

In this paper, we have discussed service by both servers in stages while moving in the forward direction. This could be extended to the case in which the service processes are in stages and can make several forward and backward jumps. In addition, it would interesting to study how the system evolves with a larger number of servers. In future work, we propose analysing a PH/two-server interdependent service model and an MAP/two-server interdependent service model and comparing them numerically.