Multi-Server Two-Way Communication Retrial Queue Subject to Disaster and Synchronous Working Vacation

Abstract

This research analyzes a multi-server retrial queue with two types of calls: working vacation and working breakdown. The incoming call may enter the retrial queue and attempt to seize a server after a random delay if all the servers are unavailable upon arrival. In its idle time, the server makes outgoing calls. All the servers take a synchronous working vacation when the system empties after regular service. The system may fail at any time due to disasters, forcing all the calls within the service area to leave the system and causing all the main servers to fail. When the main servers fail, the repair process begins immediately. The standby servers serve arriving customers at a lower level of service during the working breakdown or working vacation. For this model, we derive an explicit expression for the stationary distribution with the help of the quasi-birth-and-death process and the matrix geometric method. Further, the formulas of various system performance indices are developed. An application example is given and several numerical experiments are performed to verify the analytical results. We also perform the comparative analysis of retrial queues with/without two-way communication and two-way communication retrial queues with/without disasters. The results reveal that the proper consideration of outgoing calls to the server can reduce the average time spent in the buffer. Furthermore, a more reliable server reduces the server idle rate.

1. Introduction

In queueing theory, retrial queues have been intensive research hotspots for a long time. These retrial queues have been widely applied in analog telephone systems, call centers, computer systems, and many problems in daily life. For an overview of retrial queues, refer to the surveys [1,2,3,4]. In most studies on retrial queues, it is generally assumed that servers only serve incoming calls. However, in various service systems, especially systems with human servers, such as call centers, the servers not only need to service incoming calls but also do some private things in their free time. Researchers can read [5,6,7,8,9] to learn more about the two-way communication retrial queue approach.

In real life, for a period of time, servers may be inaccessible for various reasons. The server’s absence could mean that the server performs some other work, undergoes maintenance checks, or simply takes a break. General models in vacation queues can be found in [10,11,12,13,14]. A vacation occurs when a server is unable to serve the next customer for some reason. Allowing servers to take time off can make the queues more accurate and adaptable in situations where people are waiting in line. There could be a loss of profits if jobs accumulate during the vacation. To avoid this, servers will serve customers at different speeds. Servers give service at a reduced speed during the vacation period, rather than stopping service completely. Such a model plays an important role in some practical systems, such as network service, internet service, data transfer service, and mail service. In recent years, many authors have conducted extensive research on the retrial queues with working vacations and have achieved rich theoretical results, such as [15,16,17,18,19,20]. In practice, multi-server vacation queues, in which servers go on vacation together (called synchronous vacation) or individually (asynchronous vacations) work better than single-server queues. But, due to their complexity, analyses of these systems appear to be limited. Vijaya Laxmi and Kassahun [21] analyzed a multi-server queue with a variant working vacation. Kadi et al. [22] analyzed the synchronous multiple and single vacation policies in a multi-server queue. Bouchentouf et al. [23] investigated a finite-capacity Markovian queue with Bernoulli feedback and impatience under a synchronous multiple-vacation policy. A multi-server repairable machine system with Bernoulli feedback, customer impatience, and multiple synchronous working vacations was proposed by [24]. Ke et al. [25], which considered a multi-server retrial model with a synchronous working vacation, vacation interruption, and impatience customers under a constant retrial policy.

The authors extensively investigated retrial queues where customers are removed from the system due to a disaster or catastrophe. Modeling specialized systems such as call centers, automated teller machines, and wireless communication networks may require different types of failures. Disasters are known to be caused by negative arrivals. The presence of negative customers can cause dangers to the system. The arrival of a negative customer induces one or more positive customers to leave the system immediately. It will also cause the server to crash for a short period of time and the service channel to fail. Once a server fails, it is sent for repair, and the repair period begins immediately. The repaired server is considered as good as a new one. Negative arrivals can be seen as viruses, suppressed signals, operating errors, or system and server disasters in neural and computer communication networks. A detailed overview of negative customers is presented in [26,27,28,29,30], which discussed the retrial queues with the concept of negative arrivals.

In many real-life situations, it is common for servers to suffer lengthy and unpredictable failures while serving customers. This will cause the server to be unavailable for some time until it is repaired. Such systems with repairable servers have been extensively studied. However, in actual situations, to cope with a possible main server failure, the system may be equipped with a standby server. Kalidass and Ramanath [31] first introduced the concept of working breakdowns. When the main server is repaired, a standby server will replace the main server and serve new customers at a slower speed. Once the repair is complete, the main server is immediately back in the system and available. By allowing standby servers to provide services, the system can handle emergencies that may arise during the repair process and optimize system utilization. In addition, a working breakdown service can reduce customer complaints accruing that need to wait for the main server to be repaired and reduce the waiting cost for customers or work. Thus, working breakdown services may be a more reasonable repair strategy for unreliable retrial queues.

In most retrial queues with two-way communication, the number of servers is almost assumed to be one. But such an assumption seems unrealistic in real life. The goal of this article is to carry out an extensive analysis of a multi-server retrial queue, subject to a synchronous working vacation with two-way communication and a disaster. This model has potential applications in telephone consultation health hotline service systems. A detailed application example is provided in Section 5.

The main contributions and innovations of this paper include

• (1). Most two-way communication retrial queues in the existing literature consider a single server. This paper proposes a multi-server, two-way communication retrial queue subject to a synchronous working vacation and a disaster.

• (2). When the supplementary variable technique for deriving the closed form of the steady-state probabilities is no longer useful, we use the matrix geometric method to perform the steady-state analysis.

• (3). We derive explicit mathematical expressions for each major performance index. Through a numerical analysis, this research provides critical managerial insights for system architects or executives to evaluate system performance indices.

• (4). We conduct a multi-objective optimization study on a two-way communication retrial queue, considering a synchronous working vacation and a server failure caused by catastrophic events. This is very useful and helpful for system architects or executives to design management policies.

• (5). We compare the retrial queue without two-way communication to the retrial queue with two-way communication in terms of the average time spent in the buffer. For both discussed retrial queues, we study the impact of the outgoing generation rate on the average time spent in the buffer. By adjusting the outgoing generation rate to analyze the average time spent in the buffer, we provide actionable insights that can help system architects or executives make a more informed decision regarding the challenge of whether a server should make an outgoing call.

• (6). We also compare the server idle rate of the two-way communication retrial queue without a disaster and the two-way communication retrial queue with a disaster. For both retrial queues discussed, we study the impact of disaster rate on server idle rates. By analyzing a server’s idle rate by adjusting the disaster rate, we provide actionable insights that can help system architects or administrators make more informed decisions.

We organized this article as follows. In Section 2, we describe the retrial queue under consideration. Section 3 is devoted to the mathematical model of the steady-state analysis. The system performance indices are discussed in Section 4. Section 5 describes an application example. Section 6 contains the numerical results in the form of tables and graphs. We conclude the work in Section 7.

2. Model Description

We consider a multi-server, two-way communication retrial queue with synchronous working vacations and negative arrivals. The detailed description of the proposed retrial queue is as follows.

• Arrival and retrial process: The inbound calls obey a Poisson arrival process with rate λ. The incoming call finding that the servers fully occupied or down will either leave the system or be routed to a retrial queue and follows a constant retrial policy that attempts to seek service after an exponential distribution with rate ν. Otherwise, the call begins its service immediately. We assume that the probability of a call joining the retrial queue when the servers are in regular busy mode is b_b, the probability of a call joining the retrial queue when the servers are on vacation is b_v, and the probability of a call joining the retrial queue when the servers are under repair is b_r.

• Regular service process: The number of main servers is c. The services provided by the main servers are defined as regular services. Whenever a new incoming call or a retrial call arrives and there is at least one idle main server, the server immediately begins normal service for the call. The idle main server, after an exponentially distributed time with mean 1/α, can make an outgoing call. The time taken to serve incoming (outgoing) calls is exponentially distributed, with rate ${\mu }_1$ (${\mu }_2$).

• Removal rule and working breakdown process: When the main servers are in operation, a disaster happens at any time. The disaster obeys an exponential distribution with rate δ. When a disaster happens, all the main servers fail and all in-service calls are removed immediately from the system, but there is no impact to the orbit. The repair process begins immediately to restore its functionality after a failure event occurs. After the repair is completed, the main servers are in operation if there is one more customer in the system, including any in service. Otherwise, the main servers go on vacation at the same time. The repair time obeys an exponential distribution with rate γ. The standby servers give service at a lower rate of service for inbound calls while the main server is repaired. The working breakdown service time obeys an exponential distribution with rate ${\mu }_w$. If there are still customers in the system when the repair is completed, the standby servers will be taken out of service and the main servers will be restarted and operate at regular service rates.

• Synchronous working vacation process: When the system is found to be empty after a service or repair is completed, all the servers will take a working vacation simultaneously. The vacation time is exponentially distributed with parameter η. The servers provide a lower service rate to incoming calls arriving during this period, following an exponentially distributed rate ${\mu }_w$.

• All the stochastic processes mentioned above are assumed to be mutually independent.

3. Mathematical Model and Steady State Solution

3.1. Mathematical Model

A mathematical model of two-way communication retrial queue subject to disaster and synchronous working vacation is constructed. Let $n_1\left(t\right)$ ($n_2\left(t\right)$) denote the number of incoming (outgoing) calls under service at time t, let $L\left(t\right)$ be the size of the retrial queue at time t, and $X\left(t\right)$ represents the server state at time t. The possible values for $X\left(t\right)$ are as follows:

$X\left(t\right)=\left\{\begin{array}{l}0, \mathrm{servers} \mathrm{are} \mathrm{on} \mathrm{vacation} \mathrm{on} \mathrm{time} t,\\ 1, \mathrm{servers} \mathrm{are} \mathrm{in} \mathrm{regular} \mathrm{service} \mathrm{period} \mathrm{at} \mathrm{time} t\\ 2, \mathrm{server} \mathrm{are} \mathrm{being} \mathrm{repaired} \mathrm{at} \mathrm{time} t.\end{array}\right.,$

Obviously, $\left\{n_1\left(t\right),n_2\left(t\right),L\left(t\right),X\left(t\right)\right\}$ is a Markov process with a state space of $\mathsf{\Omega }=\left\{\left(i,0,n,k\right),i=\mathrm{0,1},\dots ,c,n\ge 0,k=\mathrm{0,2}\right\}\cup \left\{\left(i,j,n,1\right),0\le i+j\le c,n\ge 0\right\}$.

For the systems with a large or possibly infinite number of states represented in block matrix form, only the matrix geometric method can be used to derive exact solutions. If the steady-state probability for state i is known, then due to repeated behavior, the steady-state probability for state i+1 can be established. This helps to find the recursive solution of the stationary probability. According to the lexicographical sequence for the states, we have the state transition matrix with the block–tridiagonal form, as follows:

(1)$\mathbf{Q}=\left[\begin{array}{cc}\begin{array}{cc}{\mathbf{A}}^0& {\mathbf{A}}^{+}\\ {\mathbf{A}}^{-}& \mathbf{A}\end{array}& \begin{array}{cc}& \\ {\mathbf{A}}^{+}& \end{array}\\ \begin{array}{cc}& {\mathbf{A}}^{-}\\ & \end{array}& \begin{array}{cc}\mathbf{A}& {\mathbf{A}}^{+}\\ \ddots & \ddots \end{array}\end{array}\right]$

First, we define 0 to denote a matrix of appropriate dimension with all the elements equal to zero, e_i is the column-vector of dimension i, consisting of 1’s, and I_i is the identity matrix of dimension $i\times i$. For the concerned model, these matrices are represented as

${\mathbf{A}}^{+}=\left[\begin{array}{ccc}{\mathbf{A}}_v^{+}& & \\ & \begin{array}{ccc}{\mathbf{A}}_0^{+}& & \\ & \ddots & \\ & & {\mathbf{A}}_c^{+}\end{array}& \\ & & {\mathbf{A}}_r^{+}\end{array}\right],$

${\mathbf{A}}^{-}=\left[\begin{array}{ccc}{\mathbf{A}}_v^{-}& & \\ & \begin{array}{ccc}{\mathbf{A}}_0^{-}& & \\ & \ddots & \\ & & {\mathbf{A}}_c^{-}\end{array}& \\ & & {\mathbf{A}}_r^{-}\end{array}\right],$

$\mathbf{A}=\left[\begin{array}{ccc}{\mathbf{A}}_v& \begin{array}{cc}\begin{array}{cc}\eta {\mathbf{I}}_c& \end{array}& \begin{array}{cc}& \end{array}\end{array}& 0\\ \begin{array}{c}\begin{array}{c}0\\ \end{array}\\ \begin{array}{c}\\ \end{array}\end{array}& \begin{array}{cc}\begin{array}{cc}{\mathbf{A}}_0& {\mathbf{A}}_0^u\\ {\mathbf{A}}_1^l& {\mathbf{A}}_1\end{array}& \begin{array}{cc}& \\ {\mathbf{A}}_1^u& \end{array}\\ \begin{array}{cc}& \ddots \\ & \end{array}& \begin{array}{cc}\ddots & \ddots \\ {\mathbf{A}}_c^l& {\mathbf{A}}_c\end{array}\end{array}& \begin{array}{c}\begin{array}{c}{\mathbf{B}}_{c+1}\\ {\mathbf{B}}_c \end{array}\\ \begin{array}{c}⋮\\ {\mathbf{B}}_1\end{array}\end{array}\\ 0& \begin{array}{cc}\begin{array}{cc}\gamma {\mathbf{I}}_c& 0\end{array}& \begin{array}{cc}\cdots & 0\end{array}\end{array}& {\mathbf{A}}_r\end{array}\right],$

where ${\mathbf{A}}_i^{+}$, ${\mathbf{A}}_i^{-}$, ${\mathbf{A}}_i\left(i=v,r,\mathrm{0,1},\dots ,c\right)$, ${\mathbf{A}}_i^u\left(i=\mathrm{0,1},\dots ,c-1\right)$, ${\mathbf{A}}_i^l\left(i=\mathrm{1,2},\dots ,c\right)$, and ${\mathbf{B}}_i\left(i=\mathrm{1,2},\dots ,c+1\right)$ are $\left(c-i+1\right)\times \left(c-i+1\right)$, $\left(c-i+1\right)\times \left(c-i+1\right)$, $\left(c-i+1\right)\times \left(c-i+1\right)$, $\left(c-i+1\right)\times \left(c-i\right)$, $\left(c-i+1\right)\times \left(c-i+2\right),$ and $i\times \left(c+1\right)$ matrices, respectively, with entries given by

${\mathbf{A}}_i^{+}\left[k,k^{\prime }\right]=\left\{\begin{array}{c}{b}_b\lambda , k^{\prime }=k=c-i+1,\\ b_v\lambda , k^{\prime }=k=c+1, i=v,\\ b_b\lambda , k^{\prime }=k=c+1, i=r,\\ 0, \mathrm{o}.\mathrm{w}. \end{array}\right.,$

${\mathbf{A}}_i^{-}\left[k,k^{\prime }\right]=\left\{\begin{array}{l}\nu , k^{\prime }=k+1,k=\mathrm{1,2},\dots ,c-i+1, 0\le i\le c\\ \nu , k^{\prime }=k+1,k=\mathrm{1,2},\dots ,c+1, i=r,v, \\ 0,\hspace{1em}\mathrm{o}.\mathrm{w}. \end{array}\right.,$

${\mathbf{A}}_i\left[k,k^{\prime }\right]=\left\{\begin{array}{c}\lambda , k^{\prime }=k+1,k=\mathrm{0,1},2,\dots ,c-i-1,\\ k{\mu }_1, k^{\prime }=k-1,k=\mathrm{0,1},2,\dots ,c-i, \\ -{\phi }_{i,k}, k^{\prime }=k,k=\mathrm{0,1},2,\dots ,c-i, \\ 0, \mathrm{o}.\mathrm{w}. \end{array}\right.,$

${\mathbf{A}}_i^u\left[k,k^{\prime }\right]=\left\{\begin{array}{c}\left(c-i-k\right)\alpha , k^{\prime }=k,k=\mathrm{1,2},\dots ,c-i-1,\\ 0, \mathrm{o}.\mathrm{w}. \end{array}\right.,$

${\mathbf{A}}_i^l\left[k,k^{\prime }\right]=\left\{\begin{array}{c}k{\mu }_2, k^{\prime }=k,k=\mathrm{1,2},\dots ,c-i+1,\\ 0, \mathrm{o}.\mathrm{w}. \end{array}\right.,$

${\mathbf{A}}_v\left[k,k^{\prime }\right]=\left\{\begin{array}{c}\lambda , k^{\prime }=k+1,k=\mathrm{0,1},2,\dots ,c-i-1,\\ k{\mu }_w, k^{\prime }=k-1,k=\mathrm{0,1},2,\dots ,c-i, \\ -{\mathsf{\Delta }}_{v,k}, k^{\prime }=k,k=\mathrm{0,1},2,\dots ,c-i, \\ 0, \mathrm{o}.\mathrm{w}. \end{array}\right.,$

${\mathbf{A}}_r\left[k,k^{\prime }\right]=\left\{\begin{array}{c}\lambda , k^{\prime }=k+1,k=\mathrm{0,1},2,\dots ,c-i-1,\\ k{\mu }_w, k^{\prime }=k-1,k=\mathrm{0,1},2,\dots ,c-i, \\ -{\mathsf{\Delta }}_{r,k}, k^{\prime }=k,k=\mathrm{0,1},2,\dots ,c-i, \\ 0, \mathrm{o}.\mathrm{w}. \end{array}\right.,$

${\mathbf{B}}_i\left[k,k^{\prime }\right]=\left\{\begin{array}{l}\delta , k^{\prime }=1,k=\mathrm{1,2},\dots ,i\\ 0, \mathrm{o}.\mathrm{w}. \end{array}\right.,$

where

${\phi }_{i,k}=\left\{\begin{array}{c}\lambda +\left(c-i-k\right)\alpha +k{\mu }_1+i{\mu }_2+\delta +\nu , 0\le k\le c-i-1\\ b_b\lambda +\left(c-i\right){\mu }_1+i{\mu }_2+\delta , k=c-i, \end{array}\right.,$

${\mathsf{\Delta }}_{v,k}=\left\{\begin{array}{c}\lambda +k{\mu }_w+\eta +\nu , 0\le k\le c-1\\ b_v\lambda +c{\mu }_w+\eta , k=c, \end{array}\right.,$

${\mathsf{\Delta }}_{r,k}=\left\{\begin{array}{c}\lambda +k{\mu }_w+\gamma +\nu , 0\le k\le c-1\\ b_r\lambda +c{\mu }_w+\gamma , k=c. \end{array}\right.$

The matrix ${\mathbf{A}}^0$ is the same as A, but the element of ${\mathbf{A}}^0\left[c+3, c+2\right]$ is shifted to the position of ${\mathbf{A}}^0\left[c+\mathrm{3,1}\right]$, the element of ${\mathbf{A}}^0\left[2c+3, c+2\right]$ is shifted to the position of ${\mathbf{A}}^0\left[2c+\mathrm{3,1}\right],$ and the retrial rate ν is ignored.

3.2. Steady State Solution

We solve the infinitesimal matrix Q using the matrix geometric method, which is a widely used method for computing the stationary distribution of a Markov chain with a repeating structure after a certain point.

3.2.1. Ergodicity Condition

According to Neuts [32], this construction of Q indicates that $\left\{n_1\left(t\right), n_2\left(t\right), L\left(t\right), X\left(t\right)\right\}$ is a quasi-birth-and-death process. It is necessary to solve the minimal non-negative solution of the matrix quadratic equation, as shown below, to analyze this quasi-birth-and-death process:

(2)${\mathbf{R}}^2{\mathbf{A}}^{-}+\mathbf{R}\mathbf{A}+{\mathbf{A}}^{+}=0.$

It should be noted that an explicit expression for the matrix R is difficult to obtain by solving the nonlinear Equation (2). Hence, we use the following iterative procedure to approximately calculate R in our numerical study:

1. ${\mathbf{R}}_0=\mathbf{0}$ and $n=1$;

2. ${\mathbf{R}}_n=\left({\mathbf{A}}^{+}+{\mathbf{R}}_{n-1}^2{\mathbf{A}}^{-}\right){\left(-\mathbf{A}\right)}^{-1}$ for $n\ge 1$.

When $‖{\mathbf{R}}_{n+1}-{\mathbf{R}}_n‖<\mathsf{ϵ}$, we terminate the iterative procedure and return the solution of R. This iterative procedure helps us numerically calculate the stationary probabilities, that are further used to obtain various performance indices.

According to the result of Neuts [32], to prove that a quasi-birth-and-death process is a positive recurrent, it is only necessary to prove that the spectral radius of R is smaller than 1, i.e., $sp\left(\mathbf{R}\right)<1$. The matrix R fulfills $sp\left(\mathbf{R}\right)<1$ if and only if

(3)$\mathbf{y}{\mathbf{A}}^{+}{\mathbf{e}}_{\left(c+2\right)\left(c+3\right)/2}<\mathbf{y}{\mathbf{A}}^{-}{\mathbf{e}}_{\left(c+2\right)\left(c+3\right)/2},$

where y satisfies $\mathbf{y}\left({\mathbf{A}}^{+}+\mathbf{A}+{\mathbf{A}}^{-}\right)=0$ and $\mathbf{y}{\mathbf{e}}_{\left(c+2\right)\left(c+3\right)/2}=1$.

After some algebraic manipulation, the ergodicity condition is

(4)$\left(b_b\lambda +v\right)P_B<\nu -{\displaystyle \frac{\delta }{\gamma +\delta }}{\displaystyle \frac{b_r\lambda +\nu }{{\sum }_{i=0}^c{a}_i}},$

where $P_B={\sum }_{j=0}^c{y}_{c-j,j,1}$ is the blocking probability that an arriving incoming call will be blocked in the retrial queue when the servers are in busy mode, $a_i=\left(\left(i+1\right){\mu }_w{a}_{i+1}/\left(\lambda +\nu \right)\right)+\left(\gamma {\sum }_{k=i+1}^c{a}_k/\left(\lambda +\nu \right)\right)$ for $i=\mathrm{0,1},\dots ,c-1$ and $a_c=1$.

Below, a simple explicitness is obtained from Equation (4) for the special case of $c=1$. For the matrices, we have

${\mathbf{A}}^{+}=\left[\begin{array}{ccc}\begin{array}{cc}0& 0\\ 0& b_v\lambda \end{array}& \begin{array}{ccc}0& 0 & 0 \\ 0& 0 & 0 \end{array}& \begin{array}{cc}0& 0 \\ 0& 0 \end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\end{array}& \begin{array}{ccc}0& 0 & 0\\ 0& b_b\lambda & 0\\ 0& 0& b_b\lambda \end{array}& \begin{array}{cc}0& 0 \\ 0& 0 \\ 0& 0 \end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{ccc}0 & 0 & 0\\ 0 & 0 & 0\end{array}& \begin{array}{cc}0& 0 \\ 0& b_r\lambda \end{array}\end{array}\right],$

${\mathbf{A}}^{-}=\left[\begin{array}{ccc}\begin{array}{cc}0& \nu \\ 0& 0\end{array}& \begin{array}{ccc}0& 0& 0 \\ 0& 0& 0 \end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\end{array}& \begin{array}{ccc}0& \nu & 0\\ 0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{cc}0& \nu \\ 0& 0\end{array}\end{array}\right],$

$\mathbf{A}=\left[\begin{array}{ccc}\begin{array}{cc}-{\#}_1& \lambda \\ {\mu }_w& -{\#}_2\end{array}& \begin{array}{ccc}\eta & 0 & 0\\ 0 & \eta & 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\end{array}& \begin{array}{ccc}-{\#}_3& \lambda & \alpha \\ {\mu }_1& -{\#}_4& 0\\ {\mu }_2& 0& -{\#}_5\end{array}& \begin{array}{cc}\delta & 0\\ \delta & 0\\ \delta & 0\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{ccc}\gamma & 0 & 0\\ 0 & \gamma & 0\end{array}& \begin{array}{cc}-{\#}_6& \lambda \\ {\mu }_w& -{\#}_7\end{array}\end{array}\right],$

${\mathbf{A}}^{+}+\mathbf{A}+{\mathbf{A}}^{-}=\left[\begin{array}{ccc}\begin{array}{cc}-{\#}_1& \lambda +\nu \\ {\mu }_w& -{∆}_2\end{array}& \begin{array}{ccc}\eta & 0 & 0\\ 0 & \eta & 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\end{array}& \begin{array}{ccc}-{\#}_3& \lambda +\nu & \alpha \\ {\mu }_1& -{∆}_4& 0\\ {\mu }_2& 0& -{∆}_5\end{array}& \begin{array}{cc}\delta & 0\\ \delta & 0\\ \delta & 0\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{ccc}\gamma & 0 & 0\\ 0 & \gamma & 0\end{array}& \begin{array}{cc}-{\#}_6& \lambda +\nu \\ {\mu }_w& -{∆}_7\end{array}\end{array}\right],$

where ${\#}_1=\lambda +\eta +\nu$, ${\#}_2=b_v\lambda +{\mu }_w+\eta$, ${\#}_3=\lambda +\alpha +\delta +\nu$, ${\#}_4=b_b\lambda +{\mu }_1+\delta$, ${\#}_5=b_b\lambda +{\mu }_2+\delta$, ${\#}_6=\lambda +\gamma +\nu$, ${\#}_7=b_r\lambda +{\mu }_w+\gamma$, ${∆}_2={\mu }_w+\eta$, ${∆}_4={\mu }_1+\delta$, ${∆}_5={\mu }_2+\delta$ and ${∆}_7={\mu }_w+\gamma$.

Then, we have the following results:

$y_{\mathrm{0,0},0}=y_{\mathrm{1,0},0}=0,$

$y_{\mathrm{0,1},1}={\displaystyle \frac{\alpha }{{\mu }_2+\delta }}{y}_{\mathrm{0,0},1},$

$y_{\mathrm{1,0},1}={\displaystyle \frac{\delta }{\gamma +\delta }}{\displaystyle \frac{\gamma \left(\lambda +\nu \right)}{\left({\mu }_1+\delta \right)\left(\lambda +{\mu }_w+\nu +\gamma \right)}}+{\displaystyle \frac{\lambda +\nu }{{\mu }_1+\delta }}{y}_{\mathrm{0,0},1},$

$y_{\mathrm{0,0},2}={\displaystyle \frac{\delta }{\gamma +\delta }}{\displaystyle \frac{{\mu }_w+\gamma }{\lambda +{\mu }_w+\nu +\gamma }},$

$y_{\mathrm{1,0},2}={\displaystyle \frac{\delta }{\gamma +\delta }}{\displaystyle \frac{\lambda +\nu }{\lambda +{\mu }_w+\nu +\gamma }},$

where

$y_{\mathrm{0,0},1}={\displaystyle \frac{\gamma }{\gamma +\delta }}{\displaystyle \frac{\left({\mu }_2+\delta \right)\left\{{\mu }_1\left(\lambda +\nu \right)+\left({\mu }_1+\delta \right)\left({\mu }_w+\gamma \right)\right\}}{\left(\lambda +{\mu }_w+\nu +\gamma \right)\left\{\left(\lambda +{\mu }_1+\nu +\delta \right)\left({\mu }_2+\delta \right)+\alpha \left({\mu }_1+\delta \right)\right\}}}.$

Hence, the ergodicity condition is

$b_b\lambda \left({\displaystyle \frac{\lambda +\nu }{{\mu }_1+\delta }}+{\displaystyle \frac{\alpha }{{\mu }_2+\delta }}\right)<\nu +{\displaystyle \frac{\delta }{\gamma }}\vartheta ,$

where

$\vartheta ={\displaystyle \frac{\nu \left({\mu }_1+\delta \right)\left({\mu }_w+\gamma \right)-b_r\lambda \left(\lambda +\nu \right)\left({\mu }_1+\delta \right)-b_r\lambda \gamma \left(\lambda +\nu \right)}{{\mu }_1\left(\lambda +{\mu }_w+\nu +\gamma \right)+\delta \left({\mu }_r+\gamma \right)}}\left(1+{\displaystyle \frac{\lambda +\nu }{{\mu }_1+\delta }}+{\displaystyle \frac{\alpha }{{\mu }_2+\delta }}\right).$

If we set $b_b=b_r=1$ and $\delta =0$, then the ergodicity condition is reduced as follows

$\lambda \left({\displaystyle \frac{\lambda +\nu }{{\mu }_1}}+{\displaystyle \frac{\alpha }{{\mu }_2}}\right)<\nu ,$

which is consistent with Equation (14) in [33].

3.2.2. Computation of the Stationary Distribution

If $\left(b_b\lambda +v\right)P_B<\nu -{\displaystyle \frac{\delta }{\gamma +\delta }}{\displaystyle \frac{b_r\lambda +\nu }{{\sum }_{i=0}^c{a}_i}}$, let the limiting probabilities for the stationary distribution be denoted as

${\pi }_{i,j,k}\left(n\right)=\underset{t\to \infty }{\lim }Pr\left\{n_1\left(t\right)=i, n_2\left(t\right)=j, L\left(t\right)=n, X\left(t\right)=k\right\}, \left(i, j,n, k\right)\in \mathsf{\Omega }.$

Denote the stationary probability vector of the matrix Q as π. We also partitioned the vector π as $\mathsf{\pi }=\left[\mathsf{\pi }\left(0\right),\mathsf{\pi }\left(1\right),\mathsf{\pi }\left(2\right),\dots \right]$, where $\mathsf{\pi }\left(n\right)=\left[{\mathsf{\pi }}_v\left(n\right),{\mathsf{\pi }}_0\left(n\right), \dots , {\mathsf{\pi }}_c\left(n\right), {\mathsf{\pi }}_r\left(n\right)\right]$ and

${\mathsf{\pi }}_i\left(n\right)=\left\{\begin{array}{c}\left[{\pi }_{\mathrm{0,0},0}\left(n\right), {\pi }_{\mathrm{1,0},0}\left(n\right),\dots , {\pi }_{c,\mathrm{0,0}}\left(n\right)\right],i=v, \\ \left[{\pi }_{0,i,1}\left(n\right), {\pi }_{1,,i,1}\left(n\right),\dots , {\pi }_{c-i,,i,1}\left(n\right)\right],i=\mathrm{0,1},\dots ,c,\\ \left[{\pi }_{\mathrm{0,0},2}\left(n\right), {\pi }_{\mathrm{1,0},2}\left(n\right),\dots , {\pi }_{c,\mathrm{0,2}}\left(n\right)\right],i=r. \end{array}\right..$

Using matrix geometric method, we have $\mathsf{\pi }\left(n\right)=\mathsf{\pi }\left(0\right){\mathbf{R}}^n$ for $n\ge 0,$ and $\mathsf{\pi }\left(0\right)$ meets the following equations

(5)$\mathsf{\pi }\left(0\right)\left({\mathbf{A}}^0+\mathbf{R}{\mathbf{A}}^{-}\right)=0,$

(6) $\mathsf{\pi }\left(0\right){\left({\mathbf{I}}_{\left(c+2\right)\left(c+3\right)/2}-\mathbf{R}\right)}^{-1}{\mathbf{e}}_{\left(c+2\right)\left(c+3\right)/2}=1.$

Therefore, once $\mathsf{\pi }\left(0\right)$ is obtained by solving Equations (5) and (6), all other stationary probabilities are easily obtained by using $\mathsf{\pi }\left(n\right)=\mathsf{\pi }\left(0\right){\mathbf{R}}^n$ for $n\ge 0$.

4. Performance Indices

Performance indices are used to show the qualitative behavior of our model. A numerical study has been performed to find the following measure and present the parametric analysis for a decision-making goal. The expressions of the performance indices are established by using the stationary probabilities introduced in the previous section, as follows:

• Average size of the retrial queue

$\overline{O}={\sum }_{\left(i,j, k, n\right)\in \mathsf{\Omega }}n{\pi }_{i,j,k}\left(n\right);$

• Average time spent in the retrial queue

$W={\displaystyle \frac{\overline{O}}{{\lambda }_{eff}}},$

where ${\lambda }_{eff}=b_v\lambda {\sum }_{\left(i,n\right)\in \mathsf{\Omega }}{\pi }_{i,\mathrm{0,0}}\left(n\right)+b_b\lambda {\sum }_{\left(i,jn\right)\in \mathsf{\Omega }}{\pi }_{i,j,1}\left(n\right)+b_r\lambda {\sum }_{\left(i,n\right)\in \mathsf{\Omega }}{\pi }_{i,\mathrm{0,2}}\left(n\right)$;

• Probability that the servers are on vacation

$P_v={\sum }_{\left(i,n\right)\in \mathsf{\Omega }}{\pi }_{i,\mathrm{0,0}}\left(n\right);$

• Probability that the servers are in the normal busy period

$P_b={\sum }_{\left(i,jn\right)\in \mathsf{\Omega }}{\pi }_{i,j,1}\left(n\right);$

• Probability that the servers are being repaired

$P_r={\sum }_{\left(i,n\right)\in \mathsf{\Omega }}{\pi }_{i,\mathrm{0,2}}\left(n\right);$

• The average number of incoming calls under service when the servers are in the normal busy period

$E\left[S_1\right]={\sum }_{\left(i,j,n\right)\in \mathsf{\Omega }}i{\pi }_{i,j,1}\left(n\right);$

• The average number of incoming calls under service when the servers are on vacation

$E\left[S_v\right]={\sum }_{\left(i,n\right)\in \mathsf{\Omega }}i{\pi }_{i,\mathrm{0,0}}\left(n\right);$

• The average number of incoming calls under service when the servers are being repaired

$E\left[S_r\right]={\sum }_{\left(i,n\right)\in \mathsf{\Omega }}i{\pi }_{i,\mathrm{0,2}}\left(n\right);$

• Average number of outgoing calls under service when the servers are in the normal busy period

$E\left[S_2\right]={\sum }_{\left(i,j,n\right)\in \mathsf{\Omega }}j{\pi }_{i,j,1}\left(n\right);$

• The utilization of the server

$U={\displaystyle \frac{E\left[S_1\right]+E\left[S_2\right]+E\left[S_v\right]+E\left[S_r\right]}{c}}.$

5. Application Example

This work can be used to model a health helpline center. The health helpline (a centralized information system) provides services to people through telecommunications. It covers a wide range of health-related problems such as health consultations, medical advice, etc. People who need health advice can use this helpline. Calls from potential customers requesting service are considered arrivals. Arrival calls are answered by the operator, who is responsible for establishing communication between the server and the customer. The operator needs to record the order of the calls (corresponding to a “buffer”). If all the lines are busy when arriving, the arriving call has the option of leaving a message, waiting for service, and calling again after a certain amount of time (called retrial), or leaving the system (called balking). Otherwise, a chief physician (main server) or a physician assistant (standby server) serves the customer immediately. Physician assistants provide services to customers at a lower service rate only while the chief physicians are on vacation. The system may fail due to a disaster such as no network coverage or a low telephone signal status. If a disaster occurs during the customers’ consultation time, the customers being consulted will end the call. As soon as the signal strength reaches its full capacity (repair), the system provides service like new again. To avoid frequent failures, it is strongly recommended to perform regular maintenance of the service facility. The process of taking an idle server out of service for maintenance is called a vacation. If an entire service facility undergoes maintenance at the same time, we refer to this phenomenon as synchronized vacations. During maintenance, physician assistants will provide services to customers, if any. On the other hand, to minimize the chief physician’s free time, the operator will call the customer.

One scenario of the health helpline service centers is given as follows. The health helpline service, consisting of four chief physicians (main servers) (i.e., $c=4$), can be used by people who need health advice. Suppose that the number of physician assistants (standby servers) is the same as chief physicians. Calls from customers who may request service are considered to arrive at rate $\lambda =4$. Arriving calls can be served immediately if at least one chief physician/physician assistant is available and the service rate is ${\mu }_1=3$; otherwise, they have two options. First, the probability of an arriving call leaving without any attempt depends on the server states ($1-b_b=0.2$ during a regular busy period; $1-b_r=0.4$ during a repair period; and $1-b_v=0.4$ during a working vacation period). Second, it may join the buffer and must retry again, with a retrial rate of $\nu =10$. Disasters can occur at a rate of $\delta =0.1$ due to a lack of network coverage or poor telephone reception. We assume that when the signal strength is full at a rate $\gamma =1$, the system will provide service like new again. In order to avoid frequent failures, regular maintenance of service facility is performed at a rate of $\eta =0.3$. During the maintenance, the physician assistant will serve the customers at a rate of ${\mu }_w=1$. On the other hand, in order to minimize the chief physician’s free time, the operator will call the customer at a rate of $\alpha =3,$ and the service rate is ${\mu }_2=1.5$.

6. Numerical Results Discussion and Optimization Analysis

In this section, we perform performance index evaluation experiments, an optimization analysis, and parameter optimization for the system under consideration. Numerical calculations are implemented in MATLAB R2022a software.

6.1. Sensitivity Analysis

The main goal of this part is to perform a numerical examination of the efficiency of the analytical results presented previously, while showing the influence of each system parameter $\left(c, \lambda ,{\mu }_1, \alpha , \eta , \delta , {\mu }_w, \nu , b \right)$ on the system performance indices $\left(\overline{O}, P_v, P_b, P_d, U\right)$. We consider the following cases:

Case 1. $c=\left[4;6;8\right], \lambda =\left[1:0.1:5\right], \alpha =3, \nu =10, \eta =0.3, {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=1, \delta =0.1, \gamma =1, b_b=0.8, b_r=0.6, b_v=0.6$.

Case 2. $c=\left[4;6;8\right], \lambda =4, \alpha =\left[0:0.2:5\right], \nu =10, \eta =0.3, {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=1, \delta =0.1, \gamma =1, b_b=0.8, b_r=0.6, b_v=0.6$.

Case 3. $c=\left[4;6;8\right], \lambda =4, \alpha =3, \nu =10, \eta =0.3, {\mu }_1=\left[1:0.2:5\right], {\mu }_2=1.5, {\mu }_w=1, \delta =0.1, \gamma =1, b_b=0.8, b_r=0.6, b_v=0.6$.

Case 4. $c=\left[4;6;8\right], \lambda =4, \alpha =3, \nu =10, \eta =0.3, {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=\left[0:0.1:3\right], \delta =0.1, \gamma =1, b_b=0.8, b_r=0.6, b_v=0.6$.

Case 5. $c=\left[4;6;8;12\right], \lambda =4, \alpha =3, \nu =10, \eta =\left[0:0.1:3\right], {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=1, \delta =0.1, \gamma =1, b_b=0.8, b_r=0.6, b_v=0.6$.

Case 6. $c=\left[4;6;8\right], \lambda =4, \alpha =3, \nu =10, \eta =0.3, {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=1, \delta =\left[0.1:0.1:2\right], \gamma =1, b_b=0.8, b_r=0.6, b_v=0.6$.

Case 7. $c=\left[4;6;8\right], \lambda =4, \alpha =3, \nu =\left[5:0.5:20\right], \eta =0.3, {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=1, \delta =0.1, \gamma =1, b_b=0.8, b_r=0.6, b_v=0.6$.

Case 8. $c=\left[4;6;8\right], \lambda =4, \alpha =3, \nu =10, \eta =0.3, {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=1, \delta =0.1, \gamma =1,b_b=b, b_r=b, b_v=b, b=\left[0:0.05:1\right]$.

Case 9. $c=4, \lambda =\left[1:0.1:5\right], \alpha =3, \nu =10, \eta =0.3, {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=1, \delta =0.1, \gamma =1, b_b=0.8, b_r=0.6, b_v=0.6$.

Case 10. $c=4, \lambda =4, \alpha =\left[0:0.2:5\right], \nu =10, \eta =0.3, {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=1, \delta =0.1, \gamma =1, b_b=0.8, b_r=0.6, b_v=0.6$.

Case 11. $c=4, \lambda =4, \alpha =3, \nu =10, \eta =0.3, {\mu }_1=\left[1:0.2:5\right], {\mu }_2=1.5, {\mu }_w=1, \delta =0.1, \gamma =1, b_b=0.8, b_r=0.6, b_v=0.6$.

Case 12. $c=4, \lambda =4, \alpha =3, \nu =10, \eta =0.3, {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=\left[0:0.1:3\right], \delta =0.1, \gamma =1, b_b=0.8, b_r=0.6, b_v=0.6$.

Case 13. $c=4, \lambda =4, \alpha =3, \nu =10, \eta =\left[0:0.1:3\right], {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=1, \delta =0.1, \gamma =1, b_b=0.8, b_r=0.6, b_v=0.6$.

Case 14. $c=4, \lambda =4, \alpha =3, \nu =10, \eta =0.3, {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=1, \delta =\left[0.1:0.1:2\right], \gamma =1, b_b=0.8, b_r=0.6, b_v=0.6$.

Case 15. $c=4, \lambda =4, \alpha =3, \nu =\left[5:0.5:20\right], \eta =0.3, {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=1, \delta =0.1, \gamma =1, b_b=0.8, b_r=0.6, b_v=0.6$.

Case 16. $c=4, \lambda =4, \alpha =3, \nu =10, \eta =0.3, {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=1, \delta =0.1, \gamma =1, b_b=b, b_r=b, b_v=b, b=\left[0:0.05:1\right]$.

Note that the stability condition is satisfied in each case. Figure 1 and Figure 2 show the numerical results for cases 1–16, respectively. Figure 1 shows that regardless of c, an increase in the arrival rate ($\lambda$)/, the outgoing generation rate ($\alpha$)/, or the probability of joining the buffer ($b$) results in an increase in the average buffer length at the helpline service center. This is because, a higher $\lambda$/ $\alpha$/$b$ corresponds to a higher number of arrival calls, which inevitably increases the waiting time for arriving calls. In addition, the larger the value of ${\mu }_1$/${\mu }_w$/$\delta$/$\nu$, the more the average buffer length of the health helpline center will be reduced, no matter what c is. This can be explained by the fact that increasing the service rate or improving the reliability of the chief physicians will reduce the calls’ dwell time and corresponding waiting costs, thereby providing incentives for customers who may request services to enter the system. We also can find from Figure 1 that for smaller c values, the average buffer length at the health helpline center decreases as $\eta$ increases; while for larger c values, the average buffer length at the health helpline center increases slightly with the increase of $\eta$.

It is evident from Figure 2(a)–(h)-that (i) $P_b$ and $P_v$ are increasing functions of $\lambda$, $\alpha$, and $b$, while they are decreasing functions of ${\mu }_1$, ${\mu }_w$, $\eta$, $\delta$, and $\nu$; (ii) $U$ increases as $\lambda$, $\alpha$, and $b$ increase; (iii) $P_b$ is an increasing function of $\alpha$, ${\mu }_1$, $\eta$, and $\delta$, while it is a decreasing function of ${\mu }_w$; (iv) $U$ and $P_d$ seem to be unchanged by varying $\nu$; (v) $P_d$ seems to have no effect when $\lambda$ and $b$ change.

In summary, the theoretical results and numerical calculations outlined in this article can help system architects or executives identify key variables and improve system effectiveness.

6.2. Discussion and Managerial Insights

Our results have managerial implications for decision makers, concerning specific parameters that can be improved to increase the effectiveness of the proposed system. These effects can be applied to the health helpline centers mentioned in Section 5. The change in parameter values affects the average buffer length of the health helpline center and other system performance indices. The following findings are based on the numerical simulations and sensitivity analyses performed.

• To increase productivity, chief physicians initiate outgoing calls when there are no incoming calls to the system, reducing their idle time. An excessively high outgoing generation rate ($\alpha$) will only lead to track congestion, while an outgoing generation rate that is too low will cause the attending the attending physician to be idle for too long, leading to a waste of personnel.

• When the maintenance time is long ($\eta \to 0$), the entire service facility is always under maintenance, leading to track congestion. As $\eta$ increases, the probability that the chief physician is in a regular service period increases, reducing the loss of customers, thereby improving the loss of income. However, the larger the value of $\eta$, the shorter the service facility maintenance time and the average buffer length at the health helpline center. The results show that the appropriate control of the maintenance rate can effectively reduce track congestion and achieve efficient operation.

6.3. Optimization Analysis

After determining the key system parameters, we try to improve the cost–benefit ratio, including three conflicting objective functions: the small-the-better average operating cost, the small-the-better server’s idle rate, and the small-the-better average time spent in the buffer. The cost management of the service system is directly related to service quality and profit. On the other hand, the average time spent in the buffer is the most important factor in determining customers’ satisfaction with the service quality. In addition, from a management perspective, managers may wish to minimize the proportion of time a server is idle to avoid overstaffing issues, which may lead to the increased cost of an idle server. To evaluate the average operating cost, by obtaining the average buffer length at the health helpline service center and the associated probability values, the average cost function is formulated as follows.

$AC=c_0\overline{O}+c_1{P}_b+c_2{P}_v+c_3{P}_r+c_4{\mu }_1+c_5{\mu }_2+c_6{\mu }_w+c_7\eta ;$

where

$c_0\equiv \mathrm{h}\mathrm{o}\mathrm{l}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g} \mathrm{c}\mathrm{o}\mathrm{s}\mathrm{t} \mathrm{p}\mathrm{e}\mathrm{r} \mathrm{u}\mathrm{n}\mathrm{i}\mathrm{t} \mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e} \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n} \mathrm{a} \mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l} \mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{s} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{b}\mathrm{u}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{r};$

$c_1\equiv cost per unit time \mathrm{f}\mathrm{o}\mathrm{r} \mathrm{c}\mathrm{h}\mathrm{i}\mathrm{e}\mathrm{f} \mathrm{p}\mathrm{h}\mathrm{y}\mathrm{s}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{s} \mathrm{b}\mathrm{e}\mathrm{i}\mathrm{n}\mathrm{g} \mathrm{i}\mathrm{n} \mathrm{r}\mathrm{e}\mathrm{g}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{r} \mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{i}\mathrm{c}\mathrm{e} \mathrm{p}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{d};$

$c_2\equiv cost per unit time \mathrm{f}\mathrm{o}\mathrm{r} \mathrm{c}\mathrm{h}\mathrm{i}\mathrm{e}\mathrm{f} \mathrm{p}\mathrm{h}\mathrm{y}\mathrm{s}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{s} \mathrm{b}\mathrm{e}\mathrm{i}\mathrm{n}\mathrm{g} \mathrm{i}\mathrm{n} \mathrm{m}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{e} \mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{e};$

$c_3\equiv cost per unit time \mathrm{f}\mathrm{o}\mathrm{r} \mathrm{c}\mathrm{h}\mathrm{i}\mathrm{e}\mathrm{f} \mathrm{p}\mathrm{h}\mathrm{y}\mathrm{s}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{s} \mathrm{b}\mathrm{e}\mathrm{i}\mathrm{n}\mathrm{g} \mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r} \mathrm{r}\mathrm{e}\mathrm{p}\mathrm{a}\mathrm{i}\mathrm{r};$

$c_4\equiv cost \mathrm{f}\mathrm{o}\mathrm{r} \mathrm{o}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g} \mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{i}\mathrm{c}\mathrm{e} \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} \mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e} {\mu }_1;$

$c_5\equiv cost \mathrm{f}\mathrm{o}\mathrm{r} \mathrm{o}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g} \mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{i}\mathrm{c}\mathrm{e} \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} \mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e} {\mu }_2;$

$c_6\equiv cost \mathrm{f}\mathrm{o}\mathrm{r} \mathrm{o}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g} \mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{i}\mathrm{c}\mathrm{e} \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} \mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e} {\mu }_w;$

$c_7\equiv cost \mathrm{a}\mathrm{s}\mathrm{s}\mathrm{o}\mathrm{c}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d} \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} \mathrm{p}\mathrm{a}\mathrm{r}\mathrm{a}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{r} \eta \mathrm{f}\mathrm{o}\mathrm{r} \mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{i}\mathrm{c}\mathrm{e} \mathrm{f}\mathrm{a}\mathrm{c}\mathrm{i}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y} \mathrm{i}\mathrm{n} \mathrm{m}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{e} \mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{e}.$

From a management perspective, we need to minimize the server’s idle rate ($1-U$) and the average operating cost ($AC$), but, from a quality of service perspective, we also want to minimize the average time spent in the buffer, W. Therefore, our goal is to search for the optimal number of chief physicians (c), the optimal outgoing generation rate (α), and the optimal vacation rate (η), which satisfies these requirements. To this end, the multi-objective optimization problem is proposed as follows:

Min $AC$

Min $1-U$

Min $W$

Subject to $\left(b_b\lambda +v\right)P_B<\nu -{\displaystyle \frac{\delta }{\gamma +\delta }}{\displaystyle \frac{b_r\lambda +\nu }{{\sum }_{i=0}^c{a}_i}}$, $\mathsf{\alpha }\ge 0$, $\mathsf{\eta }\ge 0$, c is a positive integer.

In order to solve this problem, researchers have proposed a variety of methods, such as the multi-objective genetic algorithm [34], multi-objective particle swarm optimization [35], etc. Because multi-objective genetic algorithms do not require the user to prioritize, scale, or weigh objectives, genetic algorithms are one of the most popular metaheuristic methods for solving multi-objective optimization problems. The results are helpful for system architects or executives under the constraints of the server’s idle rate thresholds, the average time spent in the buffer thresholds, and the acceptable average operating costs.

6.3.1. Multi-Objective Genetic Algorithm

Multi-objective genetic algorithms are based on the principles of natural selection and genetics and draw on biological evolution processes to solve multi-objective optimization problems. Multi-objective genetic algorithms process populations of potential solutions, applying genetic operators such as selection, crossover, and mutation to allow the population to evolve generation after generation to obtain better solutions. Unlike single-objective optimization problems, where a clear goal guides the search for the best solution, multi-objective problems require a balance between competing goals. This equilibrium is often represented by a set of solutions called the Pareto front. The Pareto front contains non-dominated solutions, meaning that no other solution in the group is better on all objectives simultaneously, reflecting a trade-off between objectives.

The key processes of multi-objective algorithm are depicted below.

• (1). Initialization

  1. (1.1). Create an initial population.

  2. (1.2). Set the crossover rate, mutation rate, and maximal generation time.

(2). Since the stopping conditions do not apply, do the following.

• (2.1). Create an offspring population by mutation and crossover.

• (2.2). Evaluate the offspring population.

• (2.3). Combine the current population and the offspring population.

• (2.4). Compute the rank and crowding distance for all the individuals in the population obtained from (2.3).

• (2.5). Generate a new population.

• (2.6). Go to (2.1).

6.3.2. Numerical Results

A numerical example is carried out for obtaining the joint optimal values $\left(c^{*}, {\alpha }^{*}, {\eta }^{*}\right)$. The following set of cost elements and system parameters is given: $c_0=20, c_1=100, c_2=50, c_3=45, c_4=8, c_5=10, c_6=6, c_7=5, \lambda =4, \nu =10, {\mu }_1=3, {\mu }_2=1.5, {\mu }_w=1, \delta =0.1, \gamma =1, b_b=0.8, b_r=0.2, b_v=0.6$. Figure 3 depicts the Pareto optimal set. These results provide useful guidelines for system architects and executives (call managers). Table 1 summarizes the Pareto optimal solutions. After obtaining the Pareto optimal solutions, the system architect or executive can select the appropriate parameter combination from Table 1, based on his judgment and experience. For example, he chooses the appropriate combination of parameters based on the weights in his mind, with weights of 0.7, 0.15, and 0.15. Thus, the optimal solution found by setting the weights is $\left(c^{*}=4, {\alpha }^{*}=14.706, {\eta }^{*}=8.859\right)$, and the corresponding optimal objective values are $1-U=0.0936, W=3.2532, AC=378.2686$.

In Figure 3 and Table 1, one can observe that there is an inverse relationship between $AC$ and $W$ for a smaller $1-U$. The adjusted coefficient of determination is $R^2\left(adj\right)=99.98\%$, and the p-value of the F-test of the regression equation is less than 0.01. Additionally, Figure 4 presents the histogram and the Kolmogorov–Smirnov test for the regression residuals. It shows that the regression residuals have an approximately normal distribution, which means that this regression model is suitable.

The regression equation is as follows:

$$\begin{gathered} AC/1000=-60.95-0.001\left(1/W\right)+0.065W-0.819e^{-W}-0.176\log \left(W\right)+ \\ 0.659/\left(1-U\right)+32.88\left(1-U\right)+83.55e^{-\left(1-U\right)}+10.48\log \left(1-U\right). \end{gathered}$$

In order to help understand this regression, we give a simple example below. If the maximum acceptable server’s idle rate is $1-U=0.084$ and the maximum acceptable time spent in the buffer is $W=7.068$, the minimum operating cost is approximately $Ac=1000\left(-60.95-0.001\left(1/7.068\right)+0.065\times 7.068-0.819e^{-7.068}-0.176\log \left(7.068\right)+0.659/0.084+32.88\times 0.084+83.55e^{-0.084}+10.48\log \left(0.084\right)\right)=632.5548$. The actual value of the minimum operating cost is 623.538, with an error of 1.43%. Therefore, the system architect or executive can use this regression model to find the minimum operating cost that satisfies the service quality and the utilization of servers.

6.4. Numerical Results Regarding Comparison Analysis

In practice, system architects or executives are often faced with the challenge of whether a chief physician should make outgoing calls. To provide useful information for the decision process, this section firstly compares the system without two-way communication and the system with two-way communication in terms of the average time spent in the buffer. In this experiment, displayed in Figure 5, we consider the following parameters: $c=4, \lambda =4,{\mu }_1=3, {\mu }_2=1.5, {\mu }_w=1, \nu =10, \eta =0.3, \delta =0.1, \gamma =1, b_b=0.8, b_r=0.8, b_v=0.8$. It can be seen from this figure that a smaller outgoing generation rate can reduce the average time spent in the buffer. However, the outgoing generation rate is larger, which will increase the average time spent in the buffer. Therefore, in practice, system architects or executives can appropriately consider whether the chief physicians make outgoing calls.

In addition to whether the chief physician should make a call, the system architect or executive is also concerned about the external environmental factors of the system, such as network coverage or phone reception quality. In response to this problem, we also compare a two-way communication retrial queue without disasters (system I) and a two-way communication retrial queue with disasters (system II) in terms of the idle rate of the chief physicians. As can be seen from Figure 6, the idle rate of the chief physicians in System I is always lower than the idle rate of the chief physicians in System II. Therefore, if the system architects or executives can improve the external environmental factors of the system, this will help reduce the idle rate of the chief physicians.

7. Conclusions

This research focused on the analysis and optimization of a retrial queue with two-way communication, disasters, and synchronous working vacations by using the matrix–geometric approach. The matrices involved and the expression for the stability condition were derived. At the end, a numerical formula for the stationary distribution was proposed, and a related optimization model was also discussed. We also performed a comparative analysis of retrial queues with/without two-way communication and two-way communication retrial queues with/without disasters. From the numerical analysis, we found that the proper consideration of outgoing calls to the server can reduce the average time spent in the buffer. Furthermore, a more reliable server reduces the server idle rate.

This research holds potential for further extension to investigate the case where the durations of ingoing and outgoing calls follow general distributions. In addition, this can be investigated under multiple types of outgoing calls. Furthermore, it seems interesting to investigate retrial queues that include more features, such as whether the arriving call may interrupt (pre-emptive) the service of an outgoing call or join (non-pre-emptive) the retrial queue. The warm-standby machine repair system also can be modeled by a multi-server retrial queue with two-way communication. In such a system, when a machine fails, the mechanical engineer makes an incoming call towards the equipment engineer for repair. If no equipment engineer is available, the mechanical engineer leaves his contact details (joins a virtual retrial orbit). The mechanical engineer in the retrial orbit will attempt to contact them again after a random period, exhibiting a retrial behavior. On the other hand, an idle equipment engineer, after some exponentially distributed period, may make an outgoing call. There are two types of outgoing calls: in the first, after an exponentially distributed idle period, the equipment engineer is able to make an outgoing call from the source to be maintained; in the second, the equipment engineer produces outgoing calls from the retrial orbit after some idle time.

Author Contributions

Conceptualization, T.-H.L., H.-Y.H. and F.-M.C.; methodology, T.-H.L. and H.-Y.H.; software, T.-H.L. and H.-Y.H.; validation, T.-H.L. and F.-M.C.; writing—original draft preparation, T.-H.L. and H.-Y.H.; writing—review and editing, T.-H.L. and F.-M.C.; supervision, F.-M.C. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data are available upon request from the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Table 1

Optimal solutions for the multi-objective optimization model.

${\mathit{c}}^{\mathit{*}}$	${\mathit{\alpha }}^{\mathit{*}}$	${\mathit{\eta }}^{\mathit{*}}$	$1-\mathit{U}$	$\mathit{W}$	$\mathit{A}\mathit{C}$
6	0.047	19.675	0.746	0.001	236.416
4	18.199	14.537	0.076	30.545	2034.349
4	17.805	13.865	0.077	20.394	1425.580
4	17.946	13.654	0.077	23.447	1606.636
6	0.836	2.793	0.542	0.010	150.625
9	2.164	3.377	0.404	0.008	156.941
4	17.197	11.467	0.080	12.014	913.798
4	17.609	12.359	0.078	16.984	1214.696
4	16.866	11.535	0.082	9.421	759.495
4	17.418	12.380	0.079	14.378	1059.350
4	18.002	13.642	0.077	24.817	1688.237
4	18.140	12.231	0.076	28.664	1910.633
5	18.172	11.212	0.095	1.787	302.619
4	17.377	11.696	0.079	13.891	1026.901
4	17.726	12.115	0.078	18.919	1328.856
4	17.848	13.828	0.077	21.264	1477.301
4	17.903	13.078	0.077	22.458	1544.755
4	16.415	12.402	0.084	7.068	623.538
4	17.777	13.243	0.078	19.858	1390.523
4	17.869	12.022	0.077	21.700	1494.297
4	17.985	12.242	0.077	24.369	1654.525
4	18.089	11.270	0.076	27.162	1816.263
9	5.609	0.133	0.294	0.043	136.802
4	18.105	11.953	0.076	27.604	1846.024
4	17.072	11.635	0.081	10.923	849.573
4	16.639	13.276	0.083	8.107	689.851
4	16.073	13.298	0.086	5.852	555.488
4	14.420	13.481	0.095	2.945	383.013
10	0.031	0.026	0.616	0.002	96.835
4	17.555	12.212	0.078	16.179	1165.922

References

1. Artalejo, J.; Gomez-Corral, A. Retrial Queueing Systems; Springer: Berlin/Heidelberg, Germany, 2008.

2. Kim, J.; Kim, B. A survey of retrial queueing systems. Ann. Oper. Res.; 2016; 247, pp. 3-36. [DOI: https://dx.doi.org/10.1007/s10479-015-2038-7]

3. Hanukov, G.; Barron, Y.; Yechiali, U. An M/G/1 queue with repeated orbit while in service. Mathematics; 2024; 12, 3722. [DOI: https://dx.doi.org/10.3390/math12233722]

4. Kalaiselvi, J.; Saravanarajan, M.C. Exploring non-Markovian dynamics: Stochastic analysis of a single-server retrial queue with recurrent clients and balking behavior under extended Bernoulli vacations. Heliyon; 2024; 10, e33919.

5. Ayyappan, G.; Udayageetha, J.; Somasundaram, B. Analysis of non-pre-emptive priority retrial queueing system with two-way communication, Bernoulli vacation, collisions, working breakdown, immediate feedback and reneging. Int. J. Math. Oper. Res.; 2020; 16, pp. 480-498. [DOI: https://dx.doi.org/10.1504/IJMOR.2020.108420]

6. Kumar, M.S.; Dadlani, A.; Kim, K. Performance analysis of an unreliable M/G/1 retrial queue with two-way communication. Oper. Res.; 2020; 20, pp. 2267-2280.

7. Ayyappan, G.; Gowthami, R. Analysis of MAP, PH2OA/PH1I, PH2O/1 retrial queue with vacation, feedback, two-way communication and impatient customers. Soft Comput.; 2021; 25, pp. 9811-9838. [DOI: https://dx.doi.org/10.1007/s00500-020-05318-4]

8. Lee, S.W.; Kim, B.; Kim, J. Analysis of the waiting time distribution in M/G/1 retrial queues with two way communication. Ann. Oper. Res.; 2022; 310, pp. 505-518. [DOI: https://dx.doi.org/10.1007/s10479-020-03717-2]

9. Sztrik, J.; Tóth, Á.; Pintér, Á.; Bács, Z. The effect of operation time of the server on the performance of finite-source retrial queues with two-way communications to the orbit. J. Math. Sci.; 2022; 267, pp. 196-204. [DOI: https://dx.doi.org/10.1007/s10958-022-06124-z]

10. Ke, J.C.; We, C.U.; Zhang, Z.G. Recent developments in vacation models: A short survey. Int. J. Oper. Res.; 2010; 7, pp. 3-8.

11. Chandrasekaran, V.M.; Indhira, K.; Saravanarajan, M.C.; Rajadurai, P. A survey on working vacation queueing models. Int. J. Pure Appl. Math.; 2016; 106, pp. 33-41.

12. Panta, A.P.; Ghimire, R.P.; Panthi, D.; Pant, S.R. A review of vacation queueing models in different framework. Ann. Pure Appl. Math.; 2021; 24, pp. 99-121. [DOI: https://dx.doi.org/10.22457/apam.v24n2a02849]

13. Barron, Y. A threshold policy in a Markov-modulated production system with server vacation: The case of continuous and batch supplies. Adv. Appl. Probab.; 2018; 50, pp. 1246-1274. [DOI: https://dx.doi.org/10.1017/apr.2018.59]

14. Mohammadi, A.; Salehi Rad, M.R. An M/G/1 queueing model with k sequential heterogeneous service steps and vacation in the transient state. Qual. Technol. Quant. Manag.; 2022; 19, pp. 633-647. [DOI: https://dx.doi.org/10.1080/16843703.2021.1981529]

15. Yang, D.Y.; Wu, C.H. Performance analysis and optimization of a retrial queue with working vacations and starting failures. Math. Comput. Model. Dyn. Syst.; 2019; 25, pp. 463-481. [DOI: https://dx.doi.org/10.1080/13873954.2019.1660378]

16. Li, T.; Zhang, L.; Gao, S. An M/G/1 retrial queue with balking customers and Bernoulli working vacation interruption. Qual. Technol. Quant. Manag.; 2019; 16, pp. 511-530. [DOI: https://dx.doi.org/10.1080/16843703.2018.1480264]

17. Do, N.H.; Van Do, T.; Melikov, A. Equilibrium customer behavior in the M/M/1 retrial queue with working vacations and a constant retrial rate. Oper. Res.; 2020; 20, pp. 627-646. [DOI: https://dx.doi.org/10.1007/s12351-017-0369-7]

18. Gupta, P.; Kumar, N. Performance analysis of retrial queueing model with working vacation, interruption, waiting server, breakdown and repair. J. Sci. Res.; 2021; 13, pp. 833-844. [DOI: https://dx.doi.org/10.3329/jsr.v13i3.52546]

19. Jain, M.; Dhibar, S.; Sanga, S.S. Markovian working vacation queue with imperfect service, balking and retrial. J. Ambient Intell. Humaniz. Comput.; 2022; 13, pp. 1907-1923. [DOI: https://dx.doi.org/10.1007/s12652-021-02954-y]

20. GnanaSekar, M.M.N.; Kandaiyan, I. Analysis of an M/G/1 retrial queue with delayed repair and feedback under working vacation policy with impatient customers. Symmetry; 2022; 14, 2024. [DOI: https://dx.doi.org/10.3390/sym14102024]

21. Vijaya Laxmi, P.; Kassahun, T.W. Analysis of variant working vacation queue with reneging under a multi-server environment. Int. J. Manag. Sci. Eng. Manag.; 2020; 15, pp. 130-137. [DOI: https://dx.doi.org/10.1080/17509653.2019.1653237]

22. Kadi, M.; Bouchentouf, A.A.; Yahiaoui, O.L. On a multiserver queueing system with customers’ impatience until the end of service under single and multiple vacation policies. Appl. Appl. Math. Int. J.; 2020; 15, pp. 740-763.

23. Bouchentouf, A.A.; Cherfaoui, M.; Boualem, M. Analysis and performance evaluation of Markovian feedback multi-server queueing model with vacation and impatience. Am. J. Math. Manag. Sci.; 2021; 40, pp. 261-282. [DOI: https://dx.doi.org/10.1080/01966324.2020.1842271]

24. Bouchentouf, A.A.; Boualem, M.; Yahiaoui, L.; Ahmad, H. A multi-station unreliable machine model with working vacation policy and customers’ impatience. Qual. Technol. Quant. Manag.; 2022; 19, pp. 766-796. [DOI: https://dx.doi.org/10.1080/16843703.2022.2054088]

25. Ke, J.C.; Chang, F.M.; Liu, T.H. Bi-objective optimization of a retrial model with synchronous working vacation interruption. Qual. Technol. Quant. Manag.; 2024; pp. 1-24. [DOI: https://dx.doi.org/10.1080/16843703.2024.2380953]

26. Malik, G.; Upadhyaya, S.; Sharma, R. A study on retrial G-queues under different scenarios: A review. Proceedings of the International Conference on Scientific and Natural Computing; Online, 5–6 February 2021; Algorithms for Intelligent Systems Singh, D.; Awasthi, A.K.; Zelinka, I.; Deep, D. Springer: Singapore, 2021.

27. Ammar, S.I.; Rajadurai, P. Performance analysis of preemptive priority retrial queueing system with disaster under working breakdown services. Symmetry; 2019; 1, 419. [DOI: https://dx.doi.org/10.3390/sym11030419]

28. Li, J.; Li, T. An M/G/1 G-queue with server breakdown, working vacations and Bernoulli vacation interruption. IAENG Int. J. Appl. Math.; 2020; 50, pp. 421-428.

29. Gao, S.; Zhang, D.; Dong, H.; Wang, X. Equilibrium balking strategies in the repairable M/M/1 G-retrial queue with complete removals. Probab. Eng. Informational Sci.; 2021; 35, pp. 138-157. [DOI: https://dx.doi.org/10.1017/S026996481900010X]

30. Lisovskaya, E.; Fedorova, E.; Salimzyanov, R.; Moiseeva, S. Resource retrial queue with two orbits and negative customers. Mathematics; 2022; 10, 321. [DOI: https://dx.doi.org/10.3390/math10030321]

31. Kalidass, K.; Ramanath, K. A queue with working breakdowns. Comput. Ind. Eng.; 2012; 63, pp. 779-783. [DOI: https://dx.doi.org/10.1016/j.cie.2012.04.018]

32. Neuts, M.F. Matrix-Geometric Solutions in Stochastic Models, Series in the Mathematical Sciences; Johns Hopkins University Press: Baltimore, MD, USA, 1981.

33. Phung-Duc, T.; Rogiest, W. Two way communication retrial queues with balanced call blending. Proceedings of the Analytical and Stochastic Modeling Techniques and Applications: 19th International Conference, ASMTA 2012; Grenoble, France, 4–6 June 2012; Proceedings 19 Springer: Berlin/Heidelberg, Germany, 2012; pp. 16-31.

34. Konak, A.; Coit, D.W.; Smith, A.E. Multi-objective optimization using genetic algorithms: A tutorial. Reliab. Eng. Syst. Saf.; 2006; 91, pp. 992-1007. [DOI: https://dx.doi.org/10.1016/j.ress.2005.11.018]

35. Zhang, X.; Zheng, X.; Cheng, R.; Qiu, J.; Jin, Y. A competitive mechanism based multi-objective particle swarm optimizer with fast convergence. Inf. Sci.; 2018; 427, pp. 63-76. [DOI: https://dx.doi.org/10.1016/j.ins.2017.10.037]