Unreliable Multi Server Retrial Queueing System with Reneging and Diverse Outgoing Services

A queue is formed when all available servers are busy, consisting of customers waiting for service. During the wait, customers can observe the status of the servers. Queues can be found in various settings such as supermarkets, ticket counters, hospitals and etc. Dutta and Choudhury [1] utilized simulation methods to analyse a M/M/1 queues’ performance.

A retrial queue is a queuing system where customers who arrive for service but are unable to get prompt service due to occupied servers are directed to an orbit. In the orbit, they cannot observe the status of the servers and must randomly re-arrive for service after some time has elapsed. In a classical retrial queue, the re-arrival times of customers are independent of other customers. For example, in a packet-switched network, packets being sent by the router are temporarily stored in a buffer while they wait for their turn to be transmitted. Kim and Kim [2] have made a comprehensive survey of retrial queueing systems, covering their various aspects and recent developments. Zhang et al. [3] investigated a retrial queue with varying service provision rates based on the server status and proposes an optimal model for the same. Moiseev et al. [4] presented an asymptotic diffusion for a retrial queue with multiple servers, where the service time followed a hyper-exponential distribution, providing a valuable tool for analysing and optimizing such systems. Fiems [5] reviewed retrial queueing models with varying retry durations, using mathematical and probabilistic techniques to understand their performance and behaviour in practical applications. Arivudainambi and Godhandaraman [6] established the steady-state likelihoods and system effectiveness metrics of a retrial queue, including the anticipated count of customers in the queueing system and the anticipated waiting time.

Reneging is a phenomenon in queuing theory where customers in a queue anticipate a delay in receiving service and subsequently decide to leave the queue without being served. For example, a customer facing an issue on mobile connectivity network may lodge a complaint to the customer care centre. Customers can await service, when all the servers become busy, they may leave the system without getting service. Ding et al. [7] have investigated M/M/1 queues with reneging. Singh et al. [8] have considered finite M/M/1/K with reneging considering fast and slow working phases. Wang and Zhang [9] have investigated the busy period of a queue structure with discouragement of customers.

In the context of service facilities, it is acknowledged that unexpected server interruptions may be encountered due to various reasons like excessive usage, improper handling, and other unforeseen reasons. This unforeseen event may result in the temporary cessation of service delivery. After an effective repair of the inoperable server, the system resumes its normal service operation. Chang and Wang [10] studied a single server retrial queue with unreliable server and set-up duration for service resumption of off server. Kuki et al. [11] discussed a retrial queue with customer conflicts and server breakdowns. Atencia and Galán-García [12] analysed the time spent within a queueing system that incorporates breakdowns and varied retrial durations. Liu et al. [13] discussed an optimization model for a retrial queueing system involving servers that are not consistently reliable due to incomplete coverage and a delay in restarting after a shutdown. Poongothai et al. [14] and Saravanan et al. [15] analysed a two heterogeneous unreliable retrial queue with servers prone to customer discouragement. Yiming and Guo [16] analysed the behaviour of M/M/1 queue with unreliable servers. Saravanan et al. [17] established M/M/1 retrial queue with server breakdown, optional service, bakling and feedback.

Queueing theory distinguishes between two types of multi-server systems based on the similarity of their service rates. Specifically, a system is known as homogeneous if each server within the system provides service at the same rate. On the other hand, a system is termed as heterogeneous if the servers within the system have different service rates. Yang et al. [18] explored a multi-server retrial queue considering likely starting failures of servers and Bernoulli feedback. Chang et al. [19] and Ke et al. [20] analysed unreliable multi-server retrial queues subject to customer impatience with feedback. Saravanan et al. [21] discussed a multi server unreliable retrial queue involving factors like balking, reneging and synchronised vacation.

Most queueing models are only subject to serving the incoming demands from the customer end. However, for more effective performance and for boosting the system revenue, outgoing service provision from the serving end is highly appreciable during idle times. This phenomenon where the server not only provides service to the incoming demands but also makes outgoing services is known as two way communication. The servers can avail different modes of outgoing service provision like providing service through calls, emails, chats, etc. Artalejo and Phung-Duc [22] analysed a M/M/1 retrial queueing system with two-way communication between customers and the server. Phung-Duc et al. [23] analysed a retrial queue with two-way communication, applicable to call centers with call blending, for both single and multi server cases. Sakurai and Phung-Duc [24] analysed M/M/1 retrial queues with two-way communication. Further under intense situations like high congestion, slow paced retrials and immediate connections to outgoing calls. Kumar et al. [25] considered a retrial queue where a single server engages incoming calls and outgoing calls only when free, and is prone to active failures.

Utilising the available idle time effectively will have a significant impact on the financial capital from revenue perspective. With this aim in view, a study on effective use of idle time is extremely important. This need motivates the study of outgoing services from servers’ end during idle period. Further, an investigation of different types of outgoing services to achieve a significant hike in a company’s revenue is due. However, in real life practice, many service systems are providing service through multiple servers. It is essential to probe the effectiveness of several types of outgoing services in multi-server retrial service systems. Moreover, systems are often subject to technical interruptions owing to reasons varying from long exposure to faulty components. Furthermore, customers can also lose their patience being exposed to long waiting in queues. This necessitates a probe into multi-server retrial queues by considering all the discussed phenomena.

However, the existing literature on multi-server retrial queues is currently limited, despite the practical need for such systems in today’s technological landscape. Additionally, the few existing multi-server queueing models that incorporate two-way communication have only considered a single type of outgoing service. Moreover, real-world scenarios often require the consideration of unreliable systems with various technical fallbacks. To address these gaps, this article presents an analysis on multiple server retrial queues with two-way communication that accounts for different modes of outgoing service and varying technical fault rates.

The paper is organized into several sections. Section 2 elucidates the practical reasoning behind the conceptual model. Section 3 introduces a QBD model, establishing stability criteria and deriving stationary probabilities for diverse server states using the MGM. Evaluation metrics of the model are accomplished in Section 4. Section 5 presents the numerical analysis results on efficiency measures. Finally, Section 6 offers conclusive remarks.

This paper presents a study of a multiple server retrial queue involving factors like reneging and two types of outgoing service, taking into account server failures. The arrival of incoming calls is modeled as a Poisson process with rate λ, and waiting calls in the orbit can retry their request after an independent exponential distribution having parameter γ. The orbital calls attempting service and encountering all busy servers may either depart from the system with probability $\bar{r}$ or retry service with probability r. The system has m identical servers, and the service duration for incoming calls follows an exponential distribution with rate μ₁.

When a server is available, it offers either a type 1 or type 2 outgoing service, following an exponentially distributed duration at rates η₁ and η₂, respectively. The outgoing service times of type 1 and type 2 services also follow exponential distribution with rate μ₂ and μ₃, respectively. During server operation, a breakdown may occur at any time, with exponential distribution rates α₁, α₂ and α₃ for incoming service, outgoing type-1 service and outgoing type-2 service, respectively.

Once a server breaks down, it is promptly sent for repair. If a server fails while attending to a customer, that customer will transition to a waiting area and attempt their service again after a random interval. During the outgoing service, if a server fails, the customers receiving service will leave the system. The repair time for a failed server for incoming service, outgoing type-1 and type-2 service follows exponential distributions with parameters β₁, β₂ and β₃ respectively. Furthermore, the model functions based on the assumption that the time between arrivals, the incoming and outgoing service times, retrial time and the incoming and outgoing breakdown and repair times being independent of one another.

3.1 Mathematical model

Let η₁(t) represent the count of service outlets currently busy handling incoming calls, η₂(t) represent the count of service outlets currently busy handling outgoing type-1 calls, η₃(t) represent the count of service outlets currently busy handling outgoing type-2 calls, η₄(t) represent the count of service outlets that have failed and are currently unavailable, and η₅(t) represent the count of orbital customers at time t. Then, the system can be described as a continuous Markov chain {η₁(t), η₂(t), η₃(t), η₄(t), η₅(t); t≥0} with state space $\text{ }\!\!\Theta\!\!\text{ }=\left\{ \left( i,p,q,l,j \right);0\le i\le m,~0\le p\le m-i,~0\le q\le m-i, \right.~\left. 0\le l\le m-\left( i+p+q \right)~\text{and}~j\ge 0 \right\}$. The system is referred to be in the state (i, p, q, l, j) at time t when η₁(t)=i, η₂(t)=p, η₃(t)=q, η₄(t)=l and η₅(t)=j.

The steady state probability of the queueing model is expressed as $P_{p,q}^{i,l}\left( j \right)=\underset{t\to \infty }{\mathop{\text{lim}}}\,~P\left\{ {{\eta }_{1}}\left( t \right) \right.\left. =i,{{\eta }_{2}}\left( t \right)=p,~{{\eta }_{3}}\left( t \right)=q,~{{\eta }_{4}}\left( t \right)=l,~{{\eta }_{5}}\left( t \right)=j \right\}$, where 0≤i+p+q+l≤m and j≥0. The model can be investigated as a QBD process considered on the state space $\text{ }\!\!\Theta\!\!\text{ }=\left\{ \left( i,p,q,l,j \right);~0\le i+p+q+l \right.\left. ~\le m,~~j\ge 0 \right\}$. The infinitesimal generator Q of the structured form:

  $Q=\left[ \begin{array}{*{35}{l}}

   {{X}_{0}} & Y & {} & {} & {} & {} & {} & {}  \\

   {{Z}_{1}} & {{X}_{1}} & Y & {} & {} & {} & {} & {}  \\

   {} & {{Z}_{2}} & {{X}_{2}} & Y & {} & {} & {} & {}  \\

   {} & {} & \ddots  & \ddots  & \ddots  & {} & {} & {}  \\

   {} & {} & {} & {{Z}_{N}} & {{X}_{N}} & Y & {} & {}  \\

   {} & {} & {} & {} & {{Z}_{N}} & {{X}_{N}} & Y & {}  \\

   {} & {} & {} & {} & {} & \ddots  & \ddots  & \ddots   \\

\end{array} \right]$                   (1)

where, the entries X_n, n≥0, Y and Z_n, n≥1 are sub-matrices of order (m+1)(m+2)(m+3)(m+4)/24 square matrices with elements.

${{X}_{n}}=\left[ \begin{array}{*{35}{l}} {{x}_{1}} & {{c}_{1}} & {{w}_{1}} & 0  \\  {{a}_{1}} & {{x}_{2}} & 0 & {{w}_{2}}  \\ {{b}_{1}} & 0 & {{x}_{3}} & {{c}_{2}}  \\ 0 & {{b}_{2}} & {{a}_{2}} & {{x}_{4}}  \\\end{array} \right],$

$Y=\left[ \begin{array}{*{35}{l}}  {{y}_{1}} & 0 & 0 & 0  \\   0 & {{y}_{2}} & 0 & 0  \\   0 & 0 & {{y}_{3}} & 0  \\   0 & 0 & 0 & {{y}_{4}}  \\\end{array} \right],$

${{Z}_{n}}=\left[ \begin{array}{*{35}{l}}   {{z}_{1}} & 0 & 0 & 0  \\   0 & {{z}_{2}} & 0 & 0  \\   0 & 0 & {{z}_{3}} & 0  \\   0 & 0 & 0 & {{z}_{4}}  \\\end{array} \right]$

x₁ is a (m+1)(m+2)/2 sub-matrix expressed as:

${{x}_{1}}=\left[ \begin{array}{*{35}{l}}  x_{0,0}^{0} & {} & {} & {} & {}  \\   e_{0,0}^{1} & x_{0,0}^{1} & {} & {} & {}  \\   {} & \ddots  & \ddots  & {} & {}  \\   {} & {} & e_{0,0}^{m-1} & x_{0,0}^{m-1} & {}  \\   {} & {} & {} & e_{0,0}^{m} & x_{0,0}^{m}  \\ \end{array} \right]$

x₂ is a m(m+1)(m+2)/6 square matrix given by:

${{x}_{2}}=\left[ \begin{array}{*{35}{l}}   x_{1}^{1} & c_{1}^{1} & {} & {} & {} & {}  \\   a_{2}^{1} & x_{2}^{1} & c_{2}^{1} & {} & {} & {}  \\   {} & a_{3}^{1} & x_{3}^{1} & c_{3}^{1} & {} & {}  \\   {} & {} & \ddots  & \ddots  & \ddots  & {}  \\   {} & {} & {} & a_{m-1}^{1} & x_{m-1}^{1} & c_{m-1}^{1}  \\   {} & {} & {} & {} & a_{m}^{1} & x_{m}^{1}  \\   {} & {} & {} & {} & {} & {}  \\ \end{array} \right]$

where,

$x_{i}^{1}=\left[ \begin{array}{*{35}{l}}   x_{i,0}^{0} & {} & {} & {} & {}  \\   e_{i,0}^{1} & x_{i,0}^{1} & {} & {} & {}  \\   {} & \ddots  & \ddots  & {} & {}  \\   {} & {} & e_{i,0}^{m-1-i} & x_{i,0}^{m-1-i} & {}  \\   {} & {} & {} & e_{i,0}^{m-i} & x_{i,0}^{m-i}  \\ \end{array} \right],$

$a_{j}^{1}=\left[ \begin{array}{*{35}{l}}  a_{j,0}^{0} & f_{j,0}^{0} & {} & {} & {}  \\   {} & a_{j,0}^{1} & f_{j,0}^{1} & {} & {}  \\   {} & {} & \ddots  & \ddots  & {}  \\   {} & {} & {} & a_{j,0}^{m-1-j} & f_{j,0}^{m-1-j}  \\\end{array} \right],$

$c_{k}^{1}=\left[ \begin{array}{*{35}{l}}   c_{k,0}^{0} & {} & {} & {}  \\   {} & c_{k,0}^{1} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & c_{k,0}^{m-1-k}  \\   {} & {} & {} & 0  \\ \end{array} \right]$,

for 1≤i≤m, 2≤j≤m and 1≤k≤m-1.

x₃ is a m(m+1)(m+2)/6 square matrix given by:

${{x}_{3}}=\left[ \begin{array}{*{35}{l}}   x_{1}^{1} & w_{1}^{1} & {} & {} & {} & {}  \\   b_{2}^{1} & x_{2}^{1} & w_{2}^{1} & {} & {} & {}  \\   {} & b_{3}^{1} & x_{3}^{1} & w_{3}^{1} & {} & {}  \\   {} & {} & \ddots  & \ddots  & \ddots  & {}  \\   {} & {} & {} & b_{m-1}^{1} & x_{m-1}^{1} & w_{m-1}^{1}  \\   {} & {} & {} & {} & b_{m}^{1} & x_{m}^{1}  \\   {} & {} & {} & {} & {} & {}  \\ \end{array} \right]$

where,

$x_{i}^{1}=\left[ \begin{array}{*{35}{l}}   x_{0,i}^{0} & {} & {} & {} & {}  \\   e_{0,i}^{1} & x_{0,i}^{1} & {} & {} & {}  \\   {} & \ddots  & \ddots  & {} & {}  \\   {} & {} & e_{0,i}^{m-1-i} & x_{0,i}^{m-1-i} & {}  \\   {} & {} & {} & e_{0,i}^{m-i} & x_{0,i}^{m-i}  \\ \end{array} \right],$

$b_{j}^{1}=\left[ \begin{array}{*{35}{l}}   b_{0,j}^{0} & g_{0,j}^{0} & {} & {} & {}  \\   {} & b_{0,j}^{1} & g_{0,j}^{1} & {} & {}  \\   {} & {} & \ddots  & \ddots  & {}  \\   {} & {} & {} & b_{0,j}^{m-1-j} & g_{0,j}^{m-1-j}  \\ \end{array} \right],$

$w_{k}^{1}=\left[ \begin{array}{*{35}{l}}   w_{0,k}^{0} & {} & {} & {}  \\   {} & w_{0,k}^{1} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & w_{0,k}^{m-1-k}  \\   {} & {} & {} & 0  \\ \end{array} \right]$,

for 1≤i≤m, 2≤j≤m and 1≤k≤m-1.

x₄ is a m(m+1)(m²+m-2)/24 square matrix given by:

${{x}_{4}}=\left[ \begin{array}{*{35}{l}}   {{A}_{1}} & {{C}_{1}} & {} & {} & {} & {}  \\   {{B}_{2}} & {{A}_{2}} & {{C}_{2}} & {} & {} & {}  \\   {} & {{B}_{2}} & {{A}_{2}} & {{C}_{2}} & {} & {}  \\   {} & {} & \ddots  & \ddots  & \ddots  & {}  \\   {} & {} & {} & {{B}_{m-2}} & {{A}_{m-2}} & {{C}_{m-2}}  \\   {} & {} & {} & {} & {{B}_{m-1}} & {{A}_{m-1}}  \\ \end{array} \right]$

where,

${{A}_{i}}=\left[ \begin{array}{*{35}{l}}   x_{1}^{i} & w_{1}^{i} & {} & {} & {}  \\   b_{2}^{i} & x_{2}^{i} & w_{2}^{i} & {} & {}  \\   {} & \ddots  & \ddots  & \ddots  & {}  \\   {} & {} & b_{m-1-i}^{i} & x_{m-1-i}^{i} & w_{m-1-i}^{i}  \\   {} & {} & {} & b_{m-i}^{i} & x_{m-i}^{i}  \\ \end{array} \right],$

${{B}_{j}}=\left[ \begin{array}{*{35}{l}}   a_{1}^{j} & {} & {} & {} & {}  \\   {} & a_{2}^{j} & {} & {} & {}  \\   {} & {} & \ddots  & {} & {}  \\   {} & {} & {} & a_{m-j}^{j} & 0  \\ \end{array} \right],$

${{C}_{k}}=\left[ \begin{array}{*{35}{l}}   c_{1}^{k} & {} & {} & {}  \\   {} & c_{2}^{k} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & c_{m-1-k}^{k}  \\   {} & {} & {} & 0  \\ \end{array} \right]$,

for 1≤i≤m-1, 2≤j≤m-1 and 1≤k≤m-2.

where,

$x_{j}^{i}=\left[ \begin{array}{*{35}{l}}   x_{i,j}^{0} & {} & {} & {} & {}  \\   e_{i,j}^{1} & x_{i,j}^{1} & {} & {} & {}  \\   {} & \ddots  & \ddots  & {} & {}  \\   {} & {} & {} & e_{i,j}^{m-i-j} & x_{i,j}^{m-i-j}  \\ \end{array} \right],$

$b_{k}^{i}=\left[ \begin{array}{*{35}{l}}   b_{i,k}^{0} & g_{i,k}^{0} & {} & {} & {}  \\   {} & b_{i,k}^{1} & g_{i,k}^{1} & {} & {}  \\   {} & {} & \ddots  & \ddots  & {}  \\   {} & {} & {} & b_{i,k}^{m-i-k} & g_{i,k}^{m-i-k}  \\ \end{array} \right],$

$w_{l}^{i}=\left[ \begin{array}{*{35}{l}}   w_{i,l}^{0} & {} & {} & {}  \\   {} & w_{i,l}^{1} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & w_{i,l}^{m-1-i-l}  \\   {} & {} & {} & 0  \\ \end{array} \right]$,

$a_{s}^{j}=\left[ \begin{array}{*{35}{l}}   a_{j,s}^{0} & f_{j,s}^{0} & {} & {} & {}  \\   {} & a_{j,s}^{1} & f_{j,s}^{1} & {} & {}  \\   {} & {} & \ddots  & \ddots  & {}  \\   {} & {} & {} & a_{j,s}^{m-j-s} & f_{j,s}^{m-j-s}  \\\end{array} \right],$

$c_{t}^{k}=\left[ \begin{array}{*{35}{l}}   c_{k,t}^{0} & {} & {} & {}  \\   {} & c_{k,t}^{1} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & c_{k,t}^{m-1-k-t}  \\ \end{array} \right]$,

for 1≤i≤m-1, 1≤j≤m-i, 2≤k≤m-i, 1≤l≤m-1-i, 1≤s≤m-j and 1≤t≤m-1-k.

The sub-matrices a₁ and b₁ are of order [m(m+1)(m+2)/6]×[(m+1)(m+2)/2]:

${{a}_{1}}=\left[ \begin{array}{*{35}{l}}   a_{1,0}^{0} & f_{1,0}^{0} & {} & {} & {}  \\   {} & a_{1,0}^{1} & f_{1,0}^{1} & {} & {}  \\   {} & {} & \ddots  & \ddots  & {}  \\   {} & {} & {} & a_{1,0}^{m-1} & f_{1,0}^{m-1}  \\ \end{array} \right],$

${{b}_{1}}=\left[ \begin{array}{*{35}{l}}   b_{0,1}^{0} & g_{0,1}^{0} & {} & {} & {}  \\   {} & b_{0,1}^{1} & g_{0,1}^{1} & {} & {}  \\   {} & {} & \ddots  & \ddots  & {}  \\   {} & {} & {} & b_{0,1}^{m-1} & g_{0,1}^{m-1}  \\ \end{array} \right]$

The sub-matrices c₁ and w₁ are of order [(m+1)(m+2)/2]×[m(m+1)(m+2)/6]:

${{c}_{1}}=\left[ \begin{array}{*{35}{l}}   c_{0,0}^{0} & {} & {} & {}  \\   {} & c_{0,0}^{1} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & c_{0,0}^{m-1}  \\ \end{array} \right],$

${{w}_{1}}=\left[ \begin{array}{*{35}{l}}   w_{0,0}^{0} & {} & {} & {}  \\   {} & w_{0,0}^{1} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & w_{0,0}^{m-1}  \\ \end{array} \right]$

The sub-matrices a₂ and b₂ are of order [m(m+1)(m²+m-2)/24]×[m(m+1)(m+2)/6]:

${{a}_{2}}=\left[ \begin{array}{*{35}{l}}   a_{1}^{1} & {} & {} & {}  \\   {} & a_{1}^{2} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & a_{1}^{m-1}  \\ \end{array} \right]$,

$a_{1}^{i}=\left[ \begin{array}{*{35}{l}}   a_{1,i}^{0} & f_{1,i}^{0} & {} & {} & {}  \\   {} & a_{1,i}^{1} & f_{1,i}^{1} & {} & {}  \\   {} & {} & \ddots  & \ddots  & {}  \\   {} & {} & {} & a_{1,i}^{m-2-i} & f_{1,i}^{m-2-i}  \\\end{array} \right];$

 $\text{ }\!\!~\!\!\text{ }1\le i\le m-1$${{b}_{2}}=\left[ \begin{array}{*{35}{l}}   b_{1}^{1} & {} & {} & {}  \\   {} & b_{2}^{1} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & b_{m-1}^{1}  \\ \end{array} \right]$,

$b_{i}^{1}=\left[ \begin{array}{*{35}{l}}   b_{i,1}^{0} & g_{i,1}^{0} & {} & {} & {}  \\   {} & b_{i,1}^{1} & g_{i,1}^{1} & {} & {}  \\   {} & {} & \ddots  & \ddots  & {}  \\   {} & {} & {} & b_{i,1}^{m-2-i} & g_{i,1}^{m-2-i}  \\\end{array} \right];\text{ }\!\!~\!\!\text{ }$

$1\le i\le m-1$

The sub-matrices c₂ and w₂ are of order[m(m+1)(m+2)/6]×[m(m+1)(m²+m-2)/24]:

${{c}_{2}}=\left[ \begin{array}{*{35}{l}}   c_{0}^{1} & {} & {} & {}  \\   {} & c_{0}^{2} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & c_{0}^{m-1}  \\\end{array} \right],c_{0}^{i}=\left[ \begin{array}{*{35}{l}}   c_{0,i}^{0} & {} & {} & {}  \\   {} & c_{0,i}^{1} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & c_{0,i}^{m-2-i}  \\\end{array} \right]$,

for 1≤i≤m-1.

 ${{w}_{2}}=\left[ \begin{array}{*{35}{l}}   w_{1}^{0} & {} & {} & {}  \\   {} & w_{2}^{0} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & w_{m-1}^{0}  \\\end{array} \right],$

$w_{i}^{0}=\left[ \begin{array}{*{35}{l}}   w_{i,0}^{0} & {} & {} & {}  \\   {} & w_{i,0}^{1} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & w_{i,0}^{m-2-i}  \\\end{array} \right];~$

$1\le i\le m-1$

where,

$a_{p,q}^{l}={{[{{a}_{ij}}]}_{\left( m+1-\left( p+q+l \right) \right)\times \left( m+2-\left( p+q+l \right) \right)}}$,

${{a}_{ij}}=\left\{ \begin{array}{*{35}{l}}   p{{\mu }_{2}}, & i=j,1\le i\le m+1-\left( p+q+l \right)  \\   0, & \text{otherwise}  \\\end{array} \right.$,

$b_{p,q}^{l}={{[{{a}_{ij}}]}_{\left( m+1-\left( p+q+l \right) \right)\times \left( m+2-\left( p+q+l \right) \right)}}$,

${{b}_{ij}}=\left\{ \begin{array}{*{35}{l}}   q{{\mu }_{3}}, & i=j,1\le i\le m+1-\left( p+q+l \right)  \\ 0, & \text{otherwise}  \\\end{array} \right.$,

$f_{p,q}^{l}={{[{{f}_{ij}}]}_{\left( m+1-\left( p+q+l \right) \right)\times \left( m+1-\left( p+q+l \right) \right)}}$,

${{f}_{ij}}=\left\{ \begin{array}{*{35}{l}}   p{{\alpha }_{2}}, & i=j,1\le i\le m+1-\left( p+q+l \right)  \\   0, & \text{otherwise}  \\ \end{array} \right.$,

$g_{p,q}^{l}={{[{{a}_{ij}}]}_{\left( m+1-\left( p+q+l \right) \right)\times \left( m+1-\left( p+q+l \right) \right)}}$,

${{g}_{ij}}=\left\{ \begin{array}{*{35}{l}}   q{{\alpha }_{3}}, & i=j,1\le i\le m+1-\left( p+q+l \right)  \\   0, & \text{otherwise}  \\\end{array} \right.$,

$c_{p,q}^{l}={{[{{a}_{ij}}]}_{\left( m+1-\left( p+q+l \right) \right)\times \left( m-\left( p+q+l \right) \right)}}$,

${{c}_{ij}}=\left\{ \begin{array}{*{35}{l}}   {{\eta }_{1}}, & i=j,1\le i\le m-\left( p+q+l \right)  \\   0, & \text{otherwise}  \\\end{array} \right.$,

$w_{p,q}^{l}={{[{{a}_{ij}}]}_{\left( m+1-\left( p+q+l \right) \right)\times \left( m-\left( p+q+l \right) \right)}}$,

${{w}_{ij}}=\left\{ \begin{array}{*{35}{l}}   {{\eta }_{2}}, & i=j,1\le i\le m-\left( p+q+l \right)  \\   0, & \text{otherwise}  \\\end{array} \right.$,

$x_{p,q}^{l}={{[{{x}_{ij}}]}_{\left( m+1-\left( p+q+l \right) \right)\times \left( m+1-\left( p+q+l \right) \right)}}$,

$\begin{aligned} & x_{i j}= \begin{cases}\lambda, & j=i+1,1 \leq i \leq m-(p+q+l) \\ j \mu_1, & i=j+1,1 \leq j \leq m-(p+q+l) \\ -\left[\lambda+\eta_1+\eta_2+(i-1)\left(\mu_1+\alpha_1\right)+p\left(\mu_2+\alpha_2\right)\right. & \\ \left. +q\left(\mu_3+\alpha_3\right)+l\left(\beta_1+\beta_2+\beta_3\right)+n \gamma r\right], & i=j, 1 \leq i \leq m-(p+q+l) \\ -\left[\lambda+(m-(p+q+l))\left(\mu_1+\alpha_1\right)+p\left(\mu_2+\alpha_2\right)\right. & \\ \left.   +q\left(\mu_3+\alpha_3\right)+l\left(\beta_1+\beta_2+\beta_3\right)+n \gamma \bar{r}\right], & i=j=m+1-(p+q+l) \\ 0, & \text { otherwise }\end{cases} \end{aligned}$,

$e_{p,q}^{l}={{[{{e}_{ij}}]}_{\left( m+1-\left( p+q+l \right) \right)\times \left( m+2-\left( p+q+l \right) \right)}}$,

${{h}_{ij}}=\left\{ \begin{array}{*{35}{l}}   l\left( {{\beta }_{1}}+{{\beta }_{2}}+{{\beta }_{3}} \right), & i=j,1\le i\le m+1-\left( p+q+l \right)  \\   0, & \text{otherwise}  \\ \end{array} \right.$.

y₁ is a (m+1)(m+2)/2 square matrix given by:

${{y}_{1}}=\left[ \begin{array}{*{35}{l}}   y_{0,0}^{0} & d_{0,0}^{0} & {} & {} & {}  \\   {} & y_{0,0}^{1} & d_{0,0}^{1} & {} & {}  \\   {} & {} & \ddots  & \ddots  & {}  \\   {} & {} & {} & y_{0,0}^{m-1} & d_{0,0}^{m-1}  \\   {} & {} & {} & {} & y_{0,0}^{m}  \\ \end{array} \right]$

y₂ is a m(m+1)(m+2)/6 square matrix given by:

${{y}_{2}}=\left[ \begin{array}{*{35}{l}}   y_{1}^{1} & {} & {} & {}  \\   {} & y_{2}^{1} & {} & {}  \\ {} & {} & \ddots  & {}  \\   {} & {} & {} & y_{m}^{1}  \\ \end{array} \right]$,

where

$y_{i}^{1}=\left[ \begin{array}{*{35}{l}}   y_{i,0}^{0} & d_{i,0}^{0} & {} & {} & {}  \\   {} & y_{i,0}^{1} & d_{i,0}^{1} & {} & {}  \\   {} & {} & \ddots  & \ddots  & {}  \\   {} & {} & {} & y_{i,0}^{m-1-i} & d_{i,0}^{m-1-i}  \\   {} & {} & {} & {} & y_{i,0}^{m-i}  \\\end{array} \right];1\le i\le m$.

y₃ is a m(m+1)(m+2)/6 square matrix given by:

${{y}_{3}}=\left[ \begin{array}{*{35}{l}}   y_{1}^{1} & {} & {} & {}  \\   {} & y_{2}^{1} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & y_{m}^{1}  \\\end{array} \right]$,

where

$y_{j}^{1}=\left[ \begin{array}{*{35}{l}}   y_{0,j}^{0} & d_{0,j}^{0} & {} & {} & {}  \\   {} & y_{0,j}^{1} & d_{0,j}^{1} & {} & {}  \\   {} & {} & \ddots  & \ddots  & {}  \\   {} & {} & {} & y_{0,j}^{m-1-j} & d_{0,j}^{m-1-j}  \\   {} & {} & {} & {} & y_{0,j}^{m-j}  \\\end{array} \right];1\le j\le m$.

y₄ is a m(m+1)(m²+m-2)/24 square matrix given by:

${{y}_{4}}=\left[ \begin{array}{*{35}{l}}   {{D}_{1}} & {} & {} & {}  \\   {} & {{D}_{2}} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & {{D}_{m-1}}  \\\end{array} \right];{{D}_{i}}=\left[ \begin{array}{*{35}{l}}   y_{1}^{i} & {} & {} & {}  \\   {} & y_{2}^{i} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & y_{m-i}^{i}  \\ \end{array} \right],$

$y_{j}^{i}=\left[ \begin{array}{*{35}{l}}   y_{i,j}^{0} & d_{i,j}^{0} & {} & {} & {}  \\   {} & y_{i,j}^{1} & d_{i,j}^{1} & {} & {}  \\   {} & {} & \ddots  & \ddots  & {}  \\ {} & {} & {} & y_{i,j}^{m-1-i-j} & d_{i,j}^{m-1-i-j}  \\   {} & {} & {} & {} & y_{i,j}^{m-i-j}  \\\end{array} \right]$

for 1≤i≤m-1 and 1≤j≤m-I, where,

$y_{p,q}^{l}={{[{{y}_{ij}}]}_{\left( m+1-\left( p+q+l \right) \right)\times \left( m+1-\left( p+q+l \right) \right)}}$,

${{y}_{ij}}=\left\{ \begin{array}{*{35}{l}}   \lambda , & i=j=m+1-\left( p+q+l \right)  \\   0, & \text{otherwise}  \\\end{array} \right.$,

$d_{p,q}^{l}={{[{{d}_{ij}}]}_{\left( m+1-\left( p+q+l \right) \right)\times \left( m-\left( p+q+l \right) \right)}}$,

${{d}_{ij}}=\left\{ \begin{array}{*{35}{l}}   j{{\alpha }_{1}}, & i=j+1,1\le j\le m-\left( p+q+l \right)  \\   0, & \text{otherwise}  \\\end{array} \right.$.

z₁ is a (m+1)(m+2)/2 square matrix given by:

${{z}_{1}}=\left[ \begin{array}{*{35}{l}}   z_{0,0}^{0} & {} & {} & {}  \\   {} & z_{0,0}^{1} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & z_{0,0}^{m}  \\ \end{array} \right]$

z₂ is a m(m+1)(m+2)/6 square matrix given by:

${{z}_{2}}=\left[ \begin{array}{*{35}{l}}   z_{1}^{1} & {} & {} & {}  \\   {} & z_{2}^{1} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & z_{m}^{1}  \\\end{array} \right]$,

where

$z_{i}^{1}=\left[ \begin{array}{*{35}{l}}   z_{i,0}^{0} & {} & {} & {} & {}  \\   {} & z_{i,0}^{1} & {} & {} & {}  \\   {} & {} & \ddots  & {} & {}  \\   {} & {} & {} & z_{i,0}^{m-1-i} & {}  \\   {} & {} & {} & {} & z_{i,0}^{m-i}  \\\end{array} \right];~1\le i\le m$.

z₃ is m(m+1)(m+2)/6 square matrix given by:

${{z}_{3}}=\left[ \begin{array}{*{35}{l}}   z_{1}^{1} & {} & {} & {}  \\   {} & z_{2}^{1} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & z_{m}^{1}  \\\end{array} \right]$,

where

$z_{j}^{1}=\left[ \begin{array}{*{35}{l}}   z_{0,j}^{0} & {} & {} & {} & {}  \\   {} & z_{0,j}^{1} & {} & {} & {}  \\   {} & {} & \ddots  & {} & {}  \\   {} & {} & {} & z_{0,j}^{m-1-j} & {}  \\   {} & {} & {} & {} & z_{0,j}^{m-j}  \\\end{array} \right];~1\le i\le m$.

z₄ is a m(m+1)(m²+m-2)/24 square matrix given by:

${{z}_{4}}=\left[ \begin{array}{*{35}{l}}   {{E}_{1}} & {} & {} & {}  \\   {} & {{E}_{2}} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & {{E}_{m-1}}  \\ \end{array} \right];{{E}_{i}}=\left[ \begin{array}{*{35}{l}}   z_{1}^{i} & {} & {} & {}  \\   {} & z_{2}^{i} & {} & {}  \\   {} & {} & \ddots  & {}  \\   {} & {} & {} & z_{m-i}^{i}  \\ \end{array} \right],$

$z_{j}^{i}=\left[ \begin{array}{*{35}{l}}   z_{i,j}^{0} & {} & {} & {} & {}  \\   {} & z_{i,j}^{1} & {} & {} & {}  \\   {} & {} & \ddots  & {} & {}  \\   {} & {} & {} & z_{i,j}^{m-1-i-j} & {}  \\ {} & {} & {} & {} & z_{j,i}^{m-i-j}  \\ \end{array} \right]$

for 1≤i≤m-1 and 1≤j≤m-I, where,

$z_{p,q}^{l}={{[{{z}_{ij}}]}_{\left( m+1-\left( p+q+l \right) \right)\times \left( m+1-\left( p+q+l \right) \right)}}$,

${{z}_{ij}}=\left\{ \begin{array}{*{35}{l}}   n\gamma r, & j=i+1,1\le i\le m-\left( p+q+l \right)  \\   n\gamma \bar{r}, & i=j=m+1-\left( p+q+l \right)  \\ 0, & \text{otherwise}  \\ \end{array} \right.$.

The ergodicity of the Markov process Q is established when the inequality xYe<xZ_Ne is satisfied. Here, 

$\begin{gathered}x=\left[x_{0,0}^{0,0}, x_{0,0}^{1,0}, \ldots, x_{0,0}^{m, 0}, x_{0,0}^{0,1}, x_{0,0}^{1,1}, \ldots, x_{0,0}^{m-1,1}, x_{0,0}^{0,2}, \ldots, x_{0,0}^{0, m}, x_{1,0}^{0,0}, x_{1,0}^{1,0}, \ldots, x_{1,0}^{m-1,0}, x_{1,0}^{0,1}, x_{1,0}^{1,1}, \ldots, x_{1,0}^{m-2,1}, x_{1,0}^{0,2}, \ldots, x_{1,0}^{0, m-1}, x_{2,0}^{0,0}, x_{2,0}^{1,0}\right. \\ \ldots, x_{2,0}^{m-2,0}, x_{2,0}^{0,1}, \ldots, x_{2,0}^{0, m-2}, x_{3,0}^{0,0}, \ldots, x_{m, 0}^{0,0}, x_{0,1}^{0,0}, x_{0,1}^{1,0}, \ldots, x_{0,1}^{m-1,0}, x_{0,1}^{0,1}, x_{0,1}^{1,1}, \ldots, x_{0,1}^{m-2,1}, x_{0,1}^{0,2}, \ldots, x_{0,1}^{0, m-1}, x_{0,2}^{0,0}, x_{0,2}^{1,0}, \ldots, x_{0,2}^{m-2,0}, x_{0,1}^{0,1}, \\ \ldots, x_{0,2}^{0, m-2}, x_{0,3}^{0,0}, \ldots, x_{0, m}^{0,0}, x_{1,1}^{0,0}, x_{1,1}^{1,0}, \ldots, x_{1,1}^{m-2,0}, x_{1,1}^{0,1}, x_{1,1}^{1,1}, \ldots, x_{1,1}^{m-3,1}, x_{1,1}^{0,2}, \ldots, x_{1,1}^{0, m-2}, x_{1,2}^{0,0}, x_{1,2}^{1,0}, \ldots, x_{1,2}^{m-3,0}, x_{1,2}^{0,1}, \ldots, x_{1,2}^{0, m-3}, \\ \left.x_{1,3}^{0,0}, \ldots, x_{1, m-1}^{0,0}, x_{2,1}^{0,0}, x_{2,1}^{1,0}, \ldots, x_{2,1}^{m-3,0}, x_{2,1}^{0,1}, \ldots, x_{2,1}^{0, m-3}, x_{2,2}^{0,0}, \ldots, x_{2, m-2}^{0,0}, x_{3,1}^{0,0}, \ldots, x_{m-1,1}^{0,0}\right]\end{gathered}$

represents the steady state distribution of F=X_N+Y+Z_N, and e=[1, 1, …, 1]^T is a vertical array of ones.

Solving the equation xF=0, we obtain:

$x_{p,q}^{i,l}=\frac{1}{i!p!q!l!\left( {{\beta }_{1}}+{{\beta }_{2}}+{{\beta }_{3}} \right)}{{\left( \frac{\lambda +N\gamma r}{{{\mu }_{1}}+{{\alpha }_{1}}} \right)}^{i}}{{\left( \frac{{{\eta }_{1}}}{{{\mu }_{2}}+{{\alpha }_{2}}} \right)}^{p}}$ $\times {{\left( \frac{{{\eta }_{2}}}{{{\mu }_{3}}+{{\alpha }_{3}}} \right)}^{q}}\left[ \frac{{{\alpha }_{1}}\left( \lambda +N\gamma r \right)}{{{\mu }_{1}}+{{\alpha }_{1}}}+\frac{{{\alpha }_{2}}{{\eta }_{1}}}{{{\mu }_{2}}+{{\alpha }_{2}}}+\frac{{{\alpha }_{3}}{{\eta }_{2}}}{{{\mu }_{3}}+{{\alpha }_{3}}} \right]x_{0,0}^{0,0}$                           (2)

where, 0≤i≤m, 0≤p≤m-i, 0≤q≤m-(i+p), 0≤l≤m-(i+p+q).

Applying the normalization condition $xe=1,x_{0,0}^{0,0}$ is obtained as:

$x_{0,0}^{0,0}=\left\{ \underset{l=0}{\overset{m}{\mathop \sum }}\,\text{}\underset{i=0}{\overset{m-l}{\mathop \sum }}\,\text{}\underset{p=1}{\overset{m-i}{\mathop \sum }}\,\text{}\underset{q=1}{\overset{m-i}{\mathop \sum }}\,\text{}\frac{1}{i!p!q!l!\left( {{\beta }_{1}}+{{\beta }_{2}}+{{\beta }_{3}} \right)} \right.$ $\times {{\left( \frac{\lambda +N\gamma r}{{{\mu }_{1}}+{{\alpha }_{1}}} \right)}^{i}}{{\left( \frac{{{\eta }_{1}}}{{{\mu }_{2}}+{{\alpha }_{2}}} \right)}^{p}}{{\left( \frac{{{\eta }_{2}}}{{{\mu }_{3}}+{{\alpha }_{3}}} \right)}^{q}}$ $\times {{\left. \left[ \frac{{{\alpha }_{1}}\left( \lambda +N\gamma r \right)}{{{\mu }_{1}}+{{\alpha }_{1}}}+\frac{{{\alpha }_{2}}{{\eta }_{1}}}{{{\mu }_{2}}+{{\alpha }_{2}}}+\frac{{{\alpha }_{3}}{{\eta }_{2}}}{{{\mu }_{3}}+{{\alpha }_{3}}} \right] \right\}}^{-1}}$                  (3)

The criteria for stability has been obtained in Neuts [26]. It has been established that the equilibrium probabilities exist only when xYe<xZ_Ne:

$\left( \lambda -N\gamma \left( \bar{r}-r \right) \right)\left[ \underset{i=0}{\overset{m-\left( p+q \right)}{\mathop \sum }}\,\text{}\underset{q=0}{\overset{m-p}{\mathop \sum }}\,\text{}\underset{p=0}{\overset{m}{\mathop \sum }}\,\text{}x_{p,q}^{i,m-\left( p+q+i \right)} \right]$    $+{{\alpha }_{1}}\underset{l=0}{\overset{m-\left( p+q+i \right)}{\mathop \sum }}\,\text{}\underset{q=0}{\overset{m-\left( p+i \right)}{\mathop \sum }}\,\text{}\underset{p=0}{\overset{m-i}{\mathop \sum }}\,\text{}\underset{i=1}{\overset{m}{\mathop \sum }}\,\text{}i.x_{p,q}^{i,l}<N\gamma r$                    (4)

3.2 Matrix-geometric approach

The stationary probability vector $\text{ }\!\!\Pi\!\!\text{ }=\left[ {{\text{ }\!\!\Pi\!\!\text{ }}_{0}},~~{{\text{ }\!\!\Pi\!\!\text{ }}_{1}},~~{{\text{ }\!\!\Pi\!\!\text{ }}_{2}},~~\ldots  \right]$ with

$\begin{gathered}\Pi=\left[\Pi_0, \Pi_1, \Pi_2, \ldots\right] \text { with } \Pi_j=\left[P_{0,0}^{0,0}(j), P_{0,0}^{1,0}(j), \ldots, P_{0,0}^{m, 0}(j), P_{0,0}^{0,1}(j), P_{0,0}^{1,1}(j), \ldots, P_{0,0}^{m-1,1}(j), P_{0,0}^{0,2}(j), \ldots, P_{0,0}^{0, m}(j), P_{1,0}^{0,0}(j), P_{1,0}^{1,0}(j), \ldots, P_{1,0}^{m-1,0}(j), P_{1,0}^{0,1}(j),\right. \\ P_{1,0}^{1,1}(j), \ldots, P_{1,0}^{m-2,1}(j), P_{1,0}^{0,2}(j), \ldots, P_{1,0}^{0, m-1}(j), P_{2,0}^{0,0}(j), P_{2,0}^{1,0}(j), \ldots, P_{2,0}^{m-2,0}(j), P_{2,0}^{0,1}(j), \ldots, P_{2,0}^{0, m-2}(j), P_{3,0}^{0,0}(j), \ldots, P_{m, 0}^{0,0}(j), \\ P_{0,1}^{0,0}(j), P_{0,1}^{1,0}(j), \ldots, P_{0,1}^{m-1,0}(j), P_{0,1}^{0,1}(j), P_{0,1}^{1,1}(j), \ldots, P_{0,1}^{m-2,1}(j), P_{0,1}^{0,2}(j), \ldots, P_{0, m}^{0,1}(j), P_{0,2}^{0,0}(j), P_{0,2}^{1,0}(j), \ldots, P_{0,2}^{m-2,0}(j), P_{0,2}^{0,1}(j), \\ \ldots, P_{0,2}^{0, m-2}(j), P_{0,3}^{0,0}(j), \ldots, P_{0, m}^{0,0}(j), P_{1,1}^{0,0}(j), P_{1,1}^{1,0}(j), \ldots, P_{1,1}^{m-2,0}(j), P_{1,1}^{0,1}(j), P_{1,1}^{1,1}(j), \ldots, P_{1,1}^{m-3,1}(j), P_{1,1}^{0,2}(j), \ldots, \\ P_{1,1}^{0, m-2}(j), P_{1,2}^{0,0}(j), P_{1,2}^{1,0}(j), \ldots, P_{1,2}^{m-3,0}(j), P_{1,2}^{0,1}(j), \ldots, P_{1,2}^{0, m-3}(j), P_{1,3}^{0,0}(j), \ldots, P_{1, m-1}^{0,0}(j), P_{2,1}^{0,0}(j), P_{2,1}^{1,0}(j), \ldots, P_{2,1}^{m-3,0}(j), P_{2,1}^{0,1}(j), \ldots, P_{2,1}^{0, m-3}(j), \\ \left.P_{2,2}^{0,0}(j), \ldots, P_{2, m-2}^{0,0}(j), P_{3,1}^{0,0}(j), \ldots, P_{m-1,1}^{0,0}(j)\right] .\end{gathered}$.

The array of equilibrium probabilities is obtained with the help of MGM:

${{\text{ }\!\!\Pi\!\!\text{ }}_{n}}={{\text{ }\!\!\Pi\!\!\text{ }}_{N}}{{R}^{n-N}},~~n\ge N$                        (5)

where, R, the rate matrix, is the minimal nonnegative solution to Y+RX_N+R²Z_N=0. An approximation of R can be obtained in light of the sequence {R_n} as given below:

R₀=0 and ${{R}_{n+1}}=-YX_{N}^{-1}-R_{n}^{2}{{Z}_{N}}X_{N}^{-1}$.  for $n\ge 0$

{R_n} follows a monotonic pattern, allowing the evaluation of R using iterative substitution, represented as$~R=\underset{n\to \infty }{\mathop{\text{lim}}}\,{{R}_{n}}.$ This approach aligns with techniques described in Neuts [26] and Latouche and Ramaswami [27].

3.3 Stationary distribution

$\text{ }\!\!\Pi\!\!\text{ }Q=0$ results in the following:

${{\text{ }\!\!\Pi\!\!\text{ }}_{0}}{{X}_{0}}+{{\text{ }\!\!\Pi\!\!\text{ }}_{1}}{{Z}_{1}}=0$                            (6)

${{\text{ }\!\!\Pi\!\!\text{ }}_{j-1}}Y+{{\text{ }\!\!\Pi\!\!\text{ }}_{j}}{{X}_{j}}+{{\text{ }\!\!\Pi\!\!\text{ }}_{j+1}}{{Z}_{j+1}}=0,~~1\le j\le N-1$                 (7)

${{\text{ }\!\!\Pi\!\!\text{ }}_{N-1}}Y+{{\text{ }\!\!\Pi\!\!\text{ }}_{N}}{{X}_{N}}+{{\text{ }\!\!\Pi\!\!\text{ }}_{N+1}}{{Z}_{N}}=0$                 (8)

${{\text{ }\!\!\Pi\!\!\text{ }}_{N}}{{R}^{j-1-N}}Y+{{\text{ }\!\!\Pi\!\!\text{ }}_{N}}{{R}^{j-N}}{{X}_{N}}+{{\text{ }\!\!\Pi\!\!\text{ }}_{N}}{{R}^{j-N+1}}{{Z}_{N}}=0$ for j≥N+1       (9)

$\underset{j=0}{\overset{\infty }{\mathop \sum }}\,\text{}{{\text{ }\!\!\Pi\!\!\text{ }}_{j}}e=1$                   (10)

By conducting algebraic operations on Eqs. (6) through (10), we derive the subsequent recursive equations:

${{\text{ }\!\!\Pi\!\!\text{ }}_{0}}={{\text{ }\!\!\Pi\!\!\text{ }}_{1}}{{Z}_{1}}{{\left( -{{X}_{0}} \right)}^{-1}}={{\text{ }\!\!\Pi\!\!\text{ }}_{1}}{{\chi }_{1}}$                   (11)

${{\text{ }\!\!\Pi\!\!\text{ }}_{j-1}}={{\text{ }\!\!\Pi\!\!\text{ }}_{j}}{{Z}_{j}}{{\left[ -\left( {{\chi }_{j-1}}Y+{{X}_{j-1}} \right) \right]}^{-1}}={{\text{ }\!\!\Pi\!\!\text{ }}_{j}}{{\chi }_{j}},$  $~2\le j\le N$                     (12)

${{\text{ }\!\!\Pi\!\!\text{ }}_{N}}\left( {{\chi }_{N}}Y+{{X}_{N}}+R{{Z}_{N}} \right)=0$                 (13)

where, ${{\chi }_{1}}={{Z}_{1}}{{\left( -{{X}_{0}} \right)}^{-1}}$ and ${{\chi }_{j}}={{Z}_{j}}{{\left[ -\left( {{\chi }_{j-1}}Y+{{X}_{j-1}} \right) \right]}^{-1}},~$2≤j≤N.

The vector ${{\text{ }\!\!\Pi\!\!\text{ }}_{N}}$ is determined by the help of above equations. The normalization condition leads to:

$\underset{j=0}{\overset{\infty }{\mathop \sum }}\,\text{}{{\text{ }\!\!\Pi\!\!\text{ }}_{j}}e={{\text{ }\!\!\Pi\!\!\text{ }}_{N}}\left[ \underset{j=1}{\overset{N}{\mathop \sum }}\,\text{}\underset{i=N}{\overset{j}{\mathop \prod }}\,\text{}{{\chi }_{i}}+{{\left( I-R \right)}^{-1}} \right]e=1$                             (14)

Upon solving Eqs. (13) and (14) through Cramer’s rule, we acquire ${{\text{ }\!\!\Pi\!\!\text{ }}_{N}}$. Additionally, the probabilities of preceding states $\left[ {{\text{ }\!\!\Pi\!\!\text{ }}_{0}},~~{{\text{ }\!\!\Pi\!\!\text{ }}_{1}},~~\ldots ,~~{{\text{ }\!\!\Pi\!\!\text{ }}_{N-1}} \right]$ are derived from Eqs. (11) and (12). Furthermore, the probabilities of subsequent states $\left[ {{\text{ }\!\!\Pi\!\!\text{ }}_{N+1}} \right.$, $\left. {{\text{ }\!\!\Pi\!\!\text{ }}_{N+2}},~~{{\text{ }\!\!\Pi\!\!\text{ }}_{N+3}},~~\ldots  \right]$.  can be determined with the help of the formula ${{\text{ }\!\!\Pi\!\!\text{ }}_{n}}={{\text{ }\!\!\Pi\!\!\text{ }}_{N}}{{R}^{n-N}},~~n\ge N+1$.