On a retrial queue with negative customers, passive breakdown, and delayed repairs

Abstract

This paper studies an M/M/1 retrial queue with negative customers, passive breakdown, and delayed repairs. Assume that the breakdown behavior of the server during idle periods is different from that during busy periods. Passive breakdowns may occur when the server is idle, due to the lack of monitoring of the server during idle periods. When the passive breakdown occurs, the server does not get repaired immediately and enters a delayed repair phase. Negative customers arrive during the busy period, which will cause the server to break down and remove the serving customers. Under steady-state conditions, we obtain explicit expressions of the probability generating functions for the steady-state distribution, together with some important performance measures for the system. In addition, we present some numerical examples to illustrate the effects of some system parameters on important performance measures and the cost function. Finally, based on the reward-cost structure, we discuss Nash equilibrium and socially optimal strategy and numerically analyze the influence of system parameters on optimal strategies and optimal social benefits.

1. Introduction

In recent years, there has been an increasing interest in modeling and analyzing retrial queues. Retrial queues have been widely applied in many real-life systems, such as telephone switching systems, mobile communication networks, random access protocols in wireless networks, and call centers. Retrial queues can reflect the customer’s service requirements, and arriving customers (users, calls, data, packets, et al.) who find the server unavailable can join a retrial group and request their services later. Jeongsim and Bara [9] studied various retrial queueing models and mainly focused on the analytical results of queue length distribution, waiting time distribution, and tail asymptotics of queue length and the waiting time distributions. Falin [4] analyzed a retrial queueing system with batch arrival and used a generating function approach to derive the distribution of the orbit length, which yields some important performance measures. Sherman and Kharoufeh [13] analyzed an unreliable M/M/1 retrial queue with infinite capacity orbit and a normal queue. Retrial customers do not rejoin the normal queue and try to access the server directly at random intervals independently of arrivals or other retrial customers until they find the server in operation and idle. Due to unpredictable factors in reality, such as limited server lifetime, external disturbances, startup failures, etc., the servers may break down during idle or busy periods. Kulkarni and Choi [11] introduced a retrial queue with breakdown and repairs and derived stability conditions and the steady-state distribution of the system. Since then, the retrial queues with breakdowns and repairs have been studied extensively from stability and reliability perspectives. Krishna et al. [10] introduced a Markovian retrial queue with two types of breakdowns in which the server breaks down at different Poisson rates during an idle or busy period. They derived some performance measures and analyzed the orbit characteristics. Based on this, Gao et al. [5] studied the M/G/1 retrial queues with two types of breakdowns and delayed repairs and obtained some important performance measures and reliability measures by the supplementary variable method.

In recent decades, a new trend has emerged in the study of queueing systems from an economic perspective, where the most fundamental problem is determining individual equilibrium and socially optimal strategies. The study of customer behavior strategy can be traced back to Naor [12], who imposed a simple linear reward-cost structure of the observable M/M/1 model. This model has subsequently been extended and supplemented by many researchers. See Hassin and Haviv [8] and Hassin [7]. Zhu et al. [16] considered equilibrium joining strategy for the almost observable case of an unreliable Mn/G/1 queue, where the arrival rate depends on the number of customers in the system. Bontali and Economou [1] considered equilibrium joining strategy for batch service queueing systems in unobservable and observable cases. Gao et al. [6] dealt with an M/M/1 retrial queue with unreliable servers from an economic point of view. Do et al. [2] studied M/M/1 retrial queues with working vacation and constant retrial rates and obtained customer equilibrium and socially optimal strategies for different information levels. Economou and Kanta [3] considered equilibrium balking strategy in a single-server repairable queueing system under two observable information levels. Zhang et al. [15] discussed equilibrium strategies in repairable M/M/1 constant retrial queues. Wang et al. [14] discussed an M/M/1 constant retrial queue with balking customers and set-up times, and they studied equilibrium strategies in the almost unobservable queue.

Motivated by the aforementioned studies, we consider an M/M/1 retrial queueing model with two different breakdowns: passive breakdown with delayed repairs and active breakdown caused by negative customers. This model has many applications. For example, this retrial queue has potential applications in packet-switched networks. If a source host wishes to send the packets to a destination host, it first sends packets to the router it is connected to and then sends packets to the destination host. If that router is available, those packets are accepted and transmitted immediately. Otherwise, due to transmission control protocol/IP network path maximum transmission unit limitations, packets are blocked by the router, in which case the blocked packets are stored in the source host’s buffer and retransmitted after some time. The router may suffer from viruses, while transferring data, causing those packets to be lost and unable to continue transmission. If it suffers from the virus breakdown, the repair takes immediately. The system breakdown may occur, while the router is idle, then the repair will delay until the next incoming packet arrives.

As another application of the retrial model, we take a telephone consultation system as an example. In such an advice scenario, the telephone operator is responsible for establishing communication between the server and the customers; the operator needs to record the call information in a registration form (corresponding to an “orbit”). When a customer calls, if the server is idle, the operator takes down the information and the advisor serves the customer immediately. On the other hand, if the server is busy, the operator tells the customer to call again after a certain period of time (called a retrial) and the customer can choose to leave a message waiting for service or leave the system. The telephone consultation system may break down. If a breakdown occurs, while a customer is being received, it will cause the customer being consulted to end the call. If the breakdown occurs, while no one is receiving a consultation, the operator cannot immediately detect the fault and leaves it in a delayed state of repair. The server continues to serve the customer after the repair is completed.

The proposed paper aims to study the stationary performance analysis and customers’ strategic behavior in the repairable retrial queues. Our contributions are as follows: under the stability condition, we construct the balance equations to obtain the steady-state probabilities of the server in different states and derive the system performance measures. Numerical examples to analyze the effects of parameters on the performance measures. We propose a cost function and find the minimum operating cost with the numerical example. Finally, we extensively analyze the customers’ equilibrium joining behavior and socially optimal strategies.

The rest of the paper is organized as follows. Section 2 describes the model in detail. Section 3 obtains the steady-state probability distribution by using the probability generating function method for steady-state analysis. Section 4 shows some important performance measures of the system. Section 5 shows the impact of system parameters on performance measures through numerical examples. And, a cost function is proposed and optimized. Section 6 analyzes the optimal strategies of the customers through the individual utility function and social benefit function. The conclusion is given in Section 7.

2. Model description

We consider an unreliable retrial queue with two types of breakdowns: passive breakdown with delayed repairs and active breakdown caused by negative customers. Customers arrive at the system according to a Poisson process with rate λ. If a customer arrives and finds the server idle, he will immediately receive the service. We assume that the customer’s service time follows an exponential distribution with parameter µ. Otherwise, arriving customers who find the server unavailable join the orbit with probability q or balk with complementary probability $1 - q$. The retrial time obeys an exponential distribution with parameter ν. The server may encounter two types of breakdowns, namely passive breakdown and active breakdown. Passive breakdown means that when the server is idle, the server breaks down according to a Poisson process with rate η. However, due to the lack of monitoring of the server during idle periods, when a passive breakdown occurs, the server is not repaired immediately and stays there until the customer arrives at the system from outside or in orbit (if available). The repair time of passive breakdown is exponentially distributed with the parameter θ, and delayed time obeys an exponential distribution with the parameter δ. After the passive breakdown occurs, the server starts to be repaired. During this repair time, the customers do not leave the system and will receive their service immediately after completing the repair. Active breakdown means that when the server is busy, the arrival of negative customers affects the system, causing the server to break down and simultaneously forcing the customer being served to leave the system. The server is repaired immediately. The negative customers arrive according to a Poisson process with rate φ. The repair time of the active breakdown obeys an exponential distribution with the parameter β.

Finally, we assume that the inter-arrival time, the service time, the breakdown time, the repair time, the delayed time, and the retrial time are independent of each other.

Let N(t) denote the number of customers in orbit at time t and I(t) denote the state of the server at time t as defined below:

\begin{equation*}I(t) = \left\{\begin{array}{l} 0,~~\text{The server is idle},\\ 1,~~\text{The server is busy},\\ 2,~~\text{The server is being repaired due to passive breakdown},\\ 3,~~\text{The server is being repaired due to negative customers},\\ 4,~~\text{The server is in delayed repair status}. \end{array} \right.\end{equation*}

Thus, the state of the system at time t can be described by the pair $\left( {N(t),I(t)} \right)$. From Markov process theory, we know that $\left\{{\left. {\left( {N(t),I(t)} \right),t \gt 0} \right\}} \right.$ constitutes a two-dimensional Markov chain with the state space $\Omega = \{(j,i),j \ge 0,i = 0,1,2,3,4\} $. The transition rate diagram of the Markov chain is shown in Figure 1.

Figure 1. Transition rate diagram of the Markov chain.

3. Steady-state analysis

In this section, we focus on the steady-state analysis of the system.

Let $\pi (j,i)$ denote the stationary joint probability at state (j, i).

\begin{equation*}\pi (j,i) = \mathop {\lim }\limits_{t \to \infty } P\left\{{N(t) = j,I(t) = i} \right\},(j,i) \in \Omega . \end{equation*}

The probability $\pi (j,i) $ satisfies the following balance equations:

(1)\begin{equation} (\lambda + \eta )\pi (0,0) = \mu \pi (0,1) + \theta \pi (0,2), \end{equation}

(2)\begin{equation} (\lambda + \eta + \nu )\pi (j,0) = \mu \pi (j,1) + \theta \pi (j,2),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \ge 1, \end{equation}

(3)\begin{equation} (\lambda q + \mu + \varphi )\pi (0,1) = \lambda \pi (0,0) + \beta \pi (0,3) + \nu \pi (1,0), \end{equation}

(4)\begin{equation} (\lambda q + \mu + \varphi )\pi (j,1) = \lambda \pi (j,0) + \beta \pi (j,3) + \nu \pi (j + 1,0) + \lambda q\pi (j - 1,1),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \ge 1, \end{equation}

(5)\begin{equation} (\lambda q + \theta )\pi (0,2) = \delta \pi (0,4), \end{equation}

(6)\begin{equation} (\lambda q + \theta )\pi (j,2) = \delta \pi (j,4) + \lambda q\pi (j - 1,2),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \ge 1, \end{equation}

(7)\begin{equation} (\lambda q + \beta )\pi (0,3) = \varphi \pi (0,1) + \varphi \pi (1,1), \end{equation}

(8)\begin{equation} (\lambda q + \beta )\pi (j,3) = \varphi \pi (j + 1,1) + \lambda q\pi (j - 1,3),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \ge 1, \end{equation}

(9)\begin{equation} (\lambda q + \delta )\pi (0,4) = \eta \pi (0,0), \end{equation}

(10)\begin{equation} (\lambda q + \delta )\pi (j,4) = \eta \pi (j,0) + \lambda q\pi (j - 1,4),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} j \ge 1. \end{equation}

In order to solve the above equations for obtaining the system steady-state distribution, we define the partial generating functions:

\begin{equation*}{G_i}(z) = \sum\limits_{j = 0}^\infty {\pi (j,i){z^j}} ,\quad |z| \lt 1,\quad i \in \{0,1,2,3,4\} .\end{equation*}

Multiplying (1)–(10) by ${z^j}$ and summing over $j {\kern 1pt} (j \ge 0$), we obtain:

(11)\begin{equation} (\lambda + \eta + \nu ){G_0}(z) = \mu {G_1}(z) + \theta {G_2}(z) + \nu \pi (0,0), \end{equation}

(12)\begin{equation} (\lambda q + \mu + \varphi ){G_1}(z) = \lambda {G_0}(z) + \beta {G_3}(z) + \lambda qz{G_1}(z) + \frac{\nu }{z}{G_0}(z) - \frac{\nu }{z}\pi (0,0), \end{equation}

(13)\begin{equation} (\lambda q + \theta ){G_2}(z) = \delta {G_4}(z) + \lambda qz{G_2}(z), \end{equation}

(14)\begin{equation} (\lambda q + \beta ){G_3}(z) = \frac{\varphi }{z}{G_1}(z) - \frac{{\varphi (1 - z)}}{z}\pi (0,1) + \lambda qz{G_3}(z), \end{equation}

(15)\begin{equation} (\lambda q + \delta ){G_4}(z) = \eta {G_0}(z) + \lambda qz{G_4}(z). \end{equation}

Theorem 3.1. The probabilities of the server being in each state are as follows:

1. (1) The probability that the server is in a normal idle state:

   (16)\begin{equation} {P_0} = {G_0}(1) = \frac{{\mu \nu \theta \delta {\kappa _1}(\lambda q + \theta )(\lambda q + \delta ) + \beta \varphi \theta \delta \mu {\kappa _2}}}{{{\kappa _2}[\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}]}}\pi (0,1). \end{equation}

2. (2) The probability that the server is busy:

   (17)\begin{equation} {P_1} = {G_1}(1) = \frac{{\beta \mu \lambda \nu {\kappa _3}(\lambda q + \theta )(\lambda q + \delta ) + \beta \varphi \theta \delta \mu {\kappa _2}(\lambda + \nu )}}{{{\kappa _2}[\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}]}}\pi (0,1). \end{equation}

3. (3) The probability that the server is in a repaired state due to a passive breakdown:

   (18)\begin{equation} {P_2} = {G_2}(1) = \frac{{\mu \nu \delta \eta {\kappa _1}(\lambda q + \theta )(\lambda q + \delta ) + \beta \varphi \delta \mu \eta {\kappa _2}}}{{{\kappa _2}[\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}]}}\pi (0,1). \end{equation}

4. (4) The probability that the server is in a repaired state due to an active breakdown:

   (19)\begin{equation} {P_3} = {G_3}(1) = \frac{{\varphi \mu \lambda \nu {\kappa _3}(\lambda q + \theta )(\lambda q + \delta ) + {\varphi ^2}\theta \delta {\kappa _2}(\lambda + \nu )}}{{{\kappa _2}[\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}]}}\pi (0,1). \end{equation}

5. (5) The probability that the server is in a delayed repair state due to delayed repair:

   (20)\begin{equation} {P_4} = {G_4}(1) = \frac{{\mu \nu \theta \eta {\kappa _1}(\lambda q + \theta )(\lambda q + \delta ) + \beta \varphi \theta \eta \mu {\kappa _2}}}{{{\kappa _2}[\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}]}}\pi (0,1), \end{equation}
   where
   \begin{equation*}{\kappa _1} = \beta (\varphi + \mu ) - \lambda q(\beta + \varphi ),\end{equation*}

   \begin{equation*}{\kappa _2} = \lambda \left[ {\lambda q(\lambda q + \eta q + \theta + \delta ) + \eta q(\theta + \delta ) + \theta \delta } \right],\end{equation*}

   \begin{equation*}{\kappa _3} = \eta q(\theta + \delta ) + \theta \delta, \end{equation*}

   \begin{equation*}{A_1} = \nu {\kappa _1}(\lambda q + \theta )(\lambda q + \delta ) + \beta \varphi {\kappa _2},\end{equation*}

   \begin{equation*}{A_2} = \lambda \mu \nu {\kappa _3}(\lambda q + \theta )(\lambda q + \delta ) + \varphi \theta \delta {\kappa _2}(\lambda + \nu ),\end{equation*}

   (21)\begin{equation} \pi (0,1) = \frac{{{\kappa _2}\left[ {\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}} \right]}}{{\mu \left[ {\theta \delta + \eta (\theta + \delta )} \right]{A_1} + (\beta + \varphi ){A_2}}}. \end{equation}

Proof. Combining from (11) and (12), we obtain

(22)\begin{equation} \left[ {\lambda (1 - z) + \eta } \right]{G_0}(z) = \left[ {(1 - z)(\mu - \lambda qz) - \varphi z} \right]{G_1}(z) + \theta {G_0}(z) + \beta z{G_3}(z). \end{equation}

Combining (13)–(15) and (22) after some algebraic operations, we have

(23)\begin{equation} {G_1}(z) = {\Theta _1} \times \frac{{\lambda q(1 - z) + \beta }}{{\beta \varphi }}{\Theta _2}{G_0}(z) + {\Theta _2}\pi (0,1), \end{equation}

where

(24)\begin{equation} {\Theta _1} = \frac{{\lambda \left[ {\lambda q(1 - z)\left[ {\lambda q(1 - z) + \eta q + \theta + \delta } \right] + {\kappa _3}} \right]}}{{\left[ {\lambda q(1 - z) + \theta } \right]\left[ {\lambda q(1 - z) + \delta } \right]}}, \end{equation}

(25)\begin{equation} {\Theta _2} = \frac{{\beta \varphi }}{{\lambda q\left[ {(1 - z)(\mu - \lambda qz) - z(\beta + \varphi )} \right] + \beta (\varphi + \mu )}}. \end{equation}

Using the above method in (15), we have

(26)\begin{equation} {G_4}(z) = \frac{\eta }{{\lambda q(1 - z) + \delta }}{G_0}(z). \end{equation}

Substituting (26) into (13), we obtain

(27)\begin{equation} {G_2}(z) = \frac{{\delta \eta }}{{\left[ {\lambda q(1 - z) + \delta } \right]\left[ {\lambda q(1 - z) + \theta } \right]}}{G_0}(z). \end{equation}

Organizing (14), we gain

(28)\begin{equation} {G_3}(z) = \frac{\varphi }{{z\left[ {\lambda q(1 - z) + \beta } \right]}}{G_1}(z) - \frac{{\varphi (1 - z)}}{{z\left[ {\lambda q(1 - z) + \beta } \right]}}\pi (0,1). \end{equation}

After substituting the above equations into (11) with some algebraic operations, we have

(29)\begin{equation} \begin{split} {G_0}(z) = \frac{\nu }{{\lambda + \eta + \nu - \mu {\Theta _1} \times \frac{{\lambda q(1 - z) + \beta }}{{\beta \varphi }}{\Theta _2} - \frac{{\theta \delta \eta }}{{\left[ {\lambda q(1 - z) + \theta } \right]\left[ {\lambda q(1 - z) + \delta } \right]}}}}\pi (0,0)\\ + \frac{{\mu {\Theta _2}}}{{\lambda + \eta + \nu - \mu {\Theta _1} \times \frac{{\lambda q(1 - z) + \beta }}{{\beta \varphi }}{\Theta _2} - \frac{{\theta \delta \eta }}{{\left[ {\lambda q(1 - z) + \theta } \right]\left[ {\lambda q(1 - z) + \delta } \right]}}}}\pi (0,1). \end{split} \end{equation}

So far, we have obtained the requested results (23)–(29), and only $\pi (0,1)$ and $\pi (0,0)$ remains to be determined, from (1), (5), and (9), we have

(30)\begin{equation} \left[ {\lambda + \eta - \frac{{\theta \delta \eta }}{{(\lambda q + \theta )(\lambda q + \delta )}}} \right]\pi (0,0) = \mu \pi (0,1). \end{equation}

Substituting z = 1 into (24)–(28), and using the normalization condition:

(31)\begin{equation} {G_0}(1) + {G_1}(1) + {G_2}(1) + {G_3}(1) + {G_4}(1) = 1, \end{equation}

we can get $\pi (0,1).$ By calculations, the above theorem can be obtained.

From (16)–(20), the system steady-state condition is

(32)\begin{equation} \theta \delta (\lambda + \nu ){\kappa _1} \gt \beta \mu \lambda {\kappa _3}. \end{equation}

4. Performance measures

In this section, we study some important performance measures of the retrial queueing system under the stability condition (32).

1. (1) The average orbit sizes ${N_i}{\kern 1pt} {\kern 1pt} (i = 0,1,2,3,4)$ when the server is idle, busy, passive breakdown, active breakdown, and delayed time, respectively, are

   (33)\begin{equation} {N_0} = {G^{\prime}_0}\left( 1 \right) = \frac{{{\omega _1}}}{{\theta \delta {\kappa _1}\left[ {\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}} \right]}}{G_0}(1) + \frac{{\lambda q\beta \varphi \mu \theta \delta {\alpha _1}}}{{{\kappa _1}\left[ {\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}} \right]}}\pi (0,1), \end{equation}

   (34)\begin{equation} \begin{split} {N_1} &= {{G^{\prime}}_1}\left( 1 \right) = \left( {\frac{{{\omega _1}(\lambda + \nu )}}{{\theta \delta \mu {\kappa _1}[\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}]}} + \frac{{\beta \mu \lambda - {\kappa _1}(\lambda + \nu )}}{{\beta {\mu ^2}}}} \right){G_0}(1)\\ &+ \left( {\frac{{\lambda q\beta \varphi \theta \delta {\alpha _1}(\lambda + \nu )}}{{\mu {\kappa _1}{\kappa _2}[\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}]}} + \frac{{\beta \varphi {\kappa _2} - \nu {\kappa _1}(\lambda q + \theta )(\lambda q + \delta )}}{{\beta \mu {\kappa _2}}}} \right)\pi (0,1), \end{split} \end{equation}

   (35)\begin{equation} \begin{split} {N_2} &= {{G^{\prime}}_2}\left( 1 \right) = \left( {\frac{{\eta {\omega _1}}}{{{\theta ^2}\delta {\kappa _1}[\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}]}} + \frac{{\lambda q\eta (\theta + \delta )}}{{{\theta ^2}\delta }}} \right){G_0}(1)\\ &+ \frac{{\lambda q\beta \varphi \mu \theta \delta {\alpha _1}}}{{{\kappa _1}[\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}]}}\pi (0,1), \end{split} \end{equation}

   (36)\begin{equation} \begin{split} {N_3} &= {{G^{\prime}}_3}(1) = \left( {\frac{{\varphi {\omega _1}(\lambda + \nu )}}{{\beta \theta \delta \mu {\kappa _1}\left[ {\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}} \right]}} + \frac{{\beta \varphi \lambda \mu - \varphi {\kappa _1}(\lambda + \nu )}}{{{\beta ^2}{\mu ^2}}}} \right){G_0}(1)\\ &+ \left( {\frac{{\mu {\omega _2} + \lambda q\beta {\varphi ^2}\theta \delta {\alpha _1}(\lambda + \nu )}}{{\mu \beta {\kappa _1}{\kappa _2}\left[ {\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}} \right]}} + \frac{{\beta \mu \varphi {\kappa _2}(\varphi + \mu ) - \varphi \mu \nu {\kappa _1}(\lambda q + \theta )(\lambda q + \delta )}}{{{\beta ^2}{\mu ^2}{\kappa _2}}}} \right)\pi (0,1), \end{split} \end{equation}

   (37)\begin{equation} {N_4} = {G^{\prime}_4}\left( 1 \right) = \left( {\frac{{\eta {\omega _1}}}{{{\delta ^2}\theta {\kappa _1}[\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}]}} + \frac{{\lambda q\eta }}{{{\delta ^2}}}} \right){G_0}(1) + \frac{{\lambda q\beta \varphi \mu \theta \delta {\alpha _1}}}{{{\kappa _1}\left( {\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}} \right)}}\pi (0,1), \end{equation}

   where ${G_0}(1)$ and $\pi (0,1)$ are given in (16) and (21), and

   \begin{equation*}{\alpha _1} = \varphi + \mu + \beta - \lambda q,\end{equation*}

   \begin{equation*}{\omega _1} = \lambda q\eta \theta \delta \kappa _1^2(\theta + \delta ) - {\lambda ^2}q\theta \delta \mu {\kappa _1}\left[ {{\kappa _3} + \beta (\theta + \delta + q\eta )} \right] + {\lambda ^2}q\beta \mu {\kappa _3}\left[ {{\kappa _1}(\theta + \delta ) + \theta \delta {\alpha _1}} \right],\end{equation*}

   \begin{equation*}{\omega _2} = \lambda \mu \nu \varphi {\kappa _3}(\lambda q + \theta )(\lambda q + \delta )(\lambda q - \beta ) + \theta \delta {\varphi ^2}{\kappa _2}(\lambda + \nu )(\lambda q - \beta ).\end{equation*}

2. (2) Let E(N) and E(L) denote the average number of customers in the orbit and in the system, respectively. Thus,

   (38)\begin{equation} E(N) = {N_0} + {N_1} + {N_2} + {N_3} + {N_4}. \end{equation}

   The average number of customers in the system E(L), is the average number of customers in orbit plus the probability that a customer is being for service. So

   (39)\begin{equation} E(L) = E(N) + {G_1}(1) + {G_2}(1) + {G_4}(1). \end{equation}

3. (3) Assuming that a tagged customer finds that the server is unavailable on arrival and decides to join the retrial orbit, his expected (conditional) waiting time in the orbit is given by

   (40)\begin{equation} \begin{split} &T(q) = \frac{{E(N)}}{{{\lambda _{ret}}}} = \frac{E(N)}{\lambda q\left( {{G_1}(1) + {G_2}(1) + {G_3}(1) + {G_4}(1)} \right)}\\ &~~~~~~= \frac{\tau _a}{\tau _b}\left( {\frac{{\eta ({\theta ^2} + {\delta ^2} + \theta \delta )}}{{{\theta ^2}{\delta ^2}}} + \frac{{(\beta + \varphi )\left( {\beta \mu \lambda - {\kappa _1}(\lambda + \nu )} \right)}}{{\lambda q{\beta ^2}{\mu ^2}}}} \right)\\ &~~~~~~+ \left( {\frac{{\eta (\theta + \delta ) + \theta \delta }}{{\theta \delta }} + \frac{{(\beta + \varphi )(\lambda + \nu )}}{{\beta \mu }}} \right)\left( {\frac{{\beta \mu \varphi \theta \delta {\alpha _1}{\kappa _2}}}{{{\kappa _1}{\tau _a}}} + \frac{{{\omega _1}}}{{\lambda q\theta \delta {\kappa _1}\left( {\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}} \right)}}} \right)\\ &~~~~~~+ \frac{1}{{\lambda q\beta \tau _a}}\left( {{\omega _2} + \frac{{{\tau _b}(\beta + \varphi )\left( {\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}} \right)}}{{\beta \theta \delta {\mu ^2}}}} \right), \end{split} \end{equation}
   where
   \begin{equation*}{\tau _a} = \mu \nu (\lambda q + \theta )(\lambda q + \delta )\left[ {\lambda {\kappa _3}(\beta + \varphi ) + \eta {\kappa _1}(\theta + \delta )} \right] + \varphi {\kappa _2}\left[ {\theta \delta (\beta + \varphi )(\lambda + \nu ) + \beta \mu \eta (\theta + \delta )} \right],\end{equation*}

   \begin{equation*}{\tau _b} = \mu \nu \theta \delta {\kappa _1}(\lambda q + \theta )(\lambda q + \delta ) + \beta \varphi \theta \delta \mu {\kappa _2}.\end{equation*}

4. (4) The busy cycle T is defined as the length of time starting when the server completes a service and the orbit is empty and ending when the server becomes idle and the orbit is empty again. Therefore, we have

   (41)\begin{equation} E(T) = \frac{{\frac{1}{\lambda }}}{{\pi (0,0)}} = \frac{{\mu \left[ {\theta \delta + \eta (\theta + \delta )} \right]{{\rm A}_1} + (\beta + \varphi ){{\rm A}_2}}}{{\lambda \mu (\lambda q + \theta )(\lambda q + \delta )\left[ {\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}} \right]}}, \end{equation}
   where ${{\rm A}_1}, {{\rm A}_2}$ are given in (21).

5. (5) In the busy cycle, the expected length of idle period, $E({T_0})$, is computed as:

   (42)\begin{equation} E({T_0}) = E(T){G_0}(1) = \frac{{\mu \nu \theta \delta {\kappa _1}(\lambda q + \theta )(\lambda q + \delta ) + \beta \varphi \theta \delta \mu {\kappa _2}}}{{\lambda \mu (\lambda q + \theta )(\lambda q + \delta )\left[ {\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}} \right]}}. \end{equation}

6. (6) In the busy cycle, the expected length of busy period, $E({T_1})$, is determined as:

   (43)\begin{equation} E({T_1}) = E(T){G_1}(1) = \frac{{\beta \mu \lambda \nu {\kappa _3}(\lambda q + \theta )(\lambda q + \delta ) + \beta \varphi \theta \delta {\kappa _2}(\lambda + \nu )}}{{\lambda \mu (\lambda q + \theta )(\lambda q + \delta )\left[ {\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}} \right]}}. \end{equation}

7. (7) In the busy cycle, the expected length of the repair period due to passive breakdown, $E({T_2})$, is calculated as:

   (44)\begin{equation} E({T_2}) = E(T){G_2}(1) = \frac{{\mu \nu \delta \eta {\kappa _1}(\lambda q + \theta )(\lambda q + \delta ) + \beta \varphi \delta \mu \eta {\kappa _2}}}{{\lambda \mu (\lambda q + \theta )(\lambda q + \delta )\left[ {\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}} \right]}}. \end{equation}

8. (8) In the busy cycle, the expected length of the repair period due to active breakdown, $E({T_3})$, is derived as:

   (45)\begin{equation} E({T_3}) = E(T){G_3}(1) = \frac{{\varphi \mu \lambda \nu {\kappa _3}(\lambda q + \theta )(\lambda q + \delta ) + {\varphi ^2}\theta \delta {\kappa _2}(\lambda + \nu )}}{{\lambda \mu (\lambda q + \theta )(\lambda q + \delta )\left[ {\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}} \right]}}. \end{equation}

9. (9) In the busy cycle, the expected length of delayed period, $E({T_4})$, is obtained as:

   (46)\begin{equation} E({T_4}) = E(T){G_4}(1) = \frac{{\mu \nu \theta \eta {\kappa _1}(\lambda q + \theta )(\lambda q + \delta ) + \beta \varphi \theta \mu \eta {\kappa _2}}}{{\lambda \mu (\lambda q + \theta )(\lambda q + \delta )\left[ {\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}} \right]}}. \end{equation}

10. (10) The probability that the server is under breakdown/repaired is given by:

    (47)\begin{align} \begin{split} {P_R} &= {G_2}(1) + {G_3}(1) + {G_4}(1)\\ &= \frac{{\mu \nu (\lambda q + \theta )(\lambda q + \delta )[\eta {\kappa _1}(\theta + \delta ) + \lambda \varphi {\kappa _3}] + \varphi {\kappa _2}[\beta \mu \eta (\theta + \delta ) + \varphi \theta \delta (\lambda + \nu )]}}{{{\kappa _2}[\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}]}}\pi (0,1). \end{split} \end{align}

    (11) The probability of the server under the operative/working state is:

    (48)\begin{equation} {P_W} = {G_0}(1) + {G_1}(1) = \frac{{\mu \nu (\lambda q + \theta )(\lambda q + \delta )(\theta \delta {\kappa _1} + \beta \lambda {\kappa _3}) + \beta \varphi \theta \delta {\kappa _2}(\lambda + \mu + \nu )}}{{{\kappa _2}[\theta \delta {\kappa _1}(\lambda + \nu ) - \beta \mu \lambda {\kappa _3}]}}\pi (0,1). \end{equation}

5. Numerical illustrations

5.1. Impact of parameters on performance measures

In this subsection, under the system stability condition (32), we show numerical examples to explain the effects of some parameters on the performance measures.

From Figure 2, we can observe that the number of customers in the system decreases with θ and increases with η. The repair time becomes short, resulting in the server getting repaired faster and can serve more customers. As the breakdown interval becomes short, the server is more likely to cause congestion in the system queue, and the average length increases.

Figure 2. Average system length E(L) versus θ for different values of $\eta {\kern 1pt} {\kern 1pt} {\kern 1pt} (q = 0.65,\lambda = 2.2,\mu = 8,\beta = 1.4,\varphi = 0.5,\delta = 3,\nu = 3)$.

Figure 3 shows that, as intuitively expected, the expected waiting time T(q) decreases as ν increases. Furthermore, for a fixed ν, T(q) decreases with δ. This is because as the waiting time for repair becomes short, the expected waiting time for customers in the system decreases.

Figure 3. Expected waiting time T(q) versus ν for different values of $\delta {\kern 1pt} {\kern 1pt} {\kern 1pt} (q = 0.6,\lambda = 3,\mu = 8,\beta = 3.5,\varphi = 1.5,\theta = 6,\eta = 2)$.

In Figure 4, the average orbit length E(N) decreases with µ and ν. The service time and retrial time become fast, which leads to a short queue length of the orbit. Figure 5 shows the busy cycle E(T) increases with λ and decreases with µ. Because the increase in arrival rates will lead to more customers joining the system, thus increasing the busy cycle. With the increase in service rate, the server serves the customers faster, leading to a decrease in the busy cycle.

Figure 4. Average orbit length E(N) versus µ for different values of $\nu~ (q = 0.4,\lambda = 4,\beta = 1.4,\varphi = 0.5,\delta = 3,\theta = 1.5,\eta = 0.1)$.

Figure 5. The busy cycle E(T) versus λ for different values of $\mu {\kern 1pt} {\kern 1pt} {\kern 1pt} (q = 0.2,\beta = 2,\varphi = 0.7,\theta = 1.8,\delta = 3,\eta = 0.3,\nu = 3)$.

Figure 6 depicts how the server breakdown/repair probability ${P_R}$ increases with φ and decreases with β. Because the breakdown rate increases, the server tends to break down more frequently, resulting in an increased probability of the server being in a breakdown/repair state. Conversely, as the repair rate increases, the server is repaired more quickly, resulting in a lower probability of the server being in a breakdown/repair state decrease.

Figure 6. Server breakdown/repair probability ${P_R}$ versus φ for different values of $\beta {\kern 1pt} {\kern 1pt} {\kern 1pt} (q = 0.6,\lambda = 4.3,\mu = 8,\theta = 1.8,\delta = 4,\eta = 0.1,\nu = 5.5)$.

Figure 7 describes the relationship between the time occupation rates ${P_i}$ with the arrival rate λ. As the arrival rate increases, the system becomes busy. Therefore, the probability of the system is idle, passive breakdown, and delayed repair decreases, while the probability of being busy and active breakdown increases.

Figure 7. The probabilities of each state ${P_i}$ versus $\lambda {\kern 1pt} {\kern 1pt} {\kern 1pt} (q = 0.4,\beta = 4.2,\varphi = 3.4,\theta = 4,\eta = 1.1,\mu = 12,\delta = 4,\nu = 8)$.

5.2. Optimization

In this section, we present an optimization analysis of the operating cost through the constructed cost function.

Define

\begin{equation*}\begin{array}{l} {C_h} =\text{Waiting cost per unit time of the customer in the system},\\ {C_s} =\text{Cost per unit time of providing services},\\ {C_e} = \text{Cost per unit time of delayed repair time},\\ {C_o} = \text{Cost per unit time of keeping the system running},\\ {C_f} = \text{Cost per unit time of server in a passive breakdown state},\\ {C_a} = \text{Cost per unit time of server being down due to negative customers},\\ {C_i} = \text{Cost per unit time of server in a delayed repair state}. \end{array}\end{equation*}

Thus, the cost function per unit time is

(49)\begin{equation} F(\delta,\mu) = {C_h}E(N) + {C_s}\mu + {C_e}\delta + {C_o}{P_1} + {C_f}{P_2} + {C_a}{P_3} + {C_i}{P_4}. \end{equation}

5.2.1. Single-objective optimization

In the following subsection, we discuss the single-objective optimization problem, which is the expected cost function per unit time illustrated in Equation (49). However, it is difficult to obtain the optimal value because the cost function is nonlinear and very complex. Therefore, we propose some numerical calculations to optimize the cost. In the following numerical experiments, we set the parameters for

\begin{equation*}{C_h} = 0.5,{C_s} = 2,{C_e} = 2.5,{C_o} = 3,{C_f} = 4,{C_a} = 4,{C_e} = 4.5.\end{equation*}

In Figure 8, we find that $F(\delta ,\mu )$ first decreases and then increases as δ and µ increase, so we can get the minimum value. As the service rate µ increases, customers’ waiting time in the system is reduced, thus reducing the cost. The constant improvement of the service rate leads to more customers joining the system, so the cost increases. The reduction in delayed repair time reduces the time it takes for the server to enter a working state, which leads to a decrease in the cost. However, as the delayed repair rate continues to increase, more customers join the system, leading to an increase in the cost.

Figure 8. The cost variation with δ and µ $ {\kern 1pt} {\kern 1pt} {\kern 1pt} (q = 0.42,\varphi = 0.5,\eta = 0.4,\lambda = 3,\beta = 2,\theta = 1.7,\nu = 1.5)$.

From Table 1, we can see that ${\delta ^*}$ and ${\mu^*}$ increase as the arrival rate λ increases, and decrease as the active repair rate β increases. Moreover, ${\delta ^*}$ decreases as the active breakdown rate φ increases, while ${\mu ^*}$ increases with respect to φ. According to the optimal cost, $F({\delta ^*},{\mu ^*})$ increases with λ and φ. The reason is that as λ increases, more customers join the system and thus the cost increases. And with the increase of φ, the frequency of breakdowns increases, leading to an increase in customers’ waiting time in the system and an increase in the cost. In addition, It can be seen that $F({\delta ^*}, {\mu ^*})$ decreases as β increases. This is because as the repair time decreases, the customers can get service faster, which leads to the decrease of the cost.

Table 1. Optimal solutions $({\delta ^*},{\mu ^*})$ and the corresponding costs $(q = 0.42,\eta = 0.4,\nu = 1.5,\theta = 1.7)$.

5.2.2. Bi-objective optimization

Most research on the optimal design of queueing models has focused on single-objective problems where a cost or profit function is the optimization objective. However, in the real world, there are many optimization problems where multiple objective functions need to be simultaneously optimized and where multiple considerations are necessary for decision-making. Expected waiting time is the most important factor for determining customer satisfaction with service quality. In queueing systems, cost is often in conflict with service quality. In this subsection, we use the non-dominated sorting genetic algorithm (NSGA-II) to find the Pareto optimal solution set satisfying both the minimum cost and the minimum waiting time. We build a bi-objective optimization model to minimize both the expected cost ${F(\delta,\mu)}$ and the expected waiting time ${T(\delta,\mu)}$, and consider the regression relationships that exist between them.

The bi-objective optimization problem is formulated as:

(50)\begin{equation} \min {\kern 1pt} {\kern 1pt} {\kern 1pt} [{F(\delta,\mu)}{\kern 1pt} {\kern 1pt} {T(\delta,\mu)}]. \end{equation}

Multi-objective genetic algorithm is an evolutionary algorithm used to analyze and solve multi-objective optimization problems that is based on the genetic algorithm and Pareto optimal concept. Its core is to coordinate the relationship between each objective function and find the optimal solution set that makes each objective function as small as possible (or relatively large). We select the following system parameters: $q = 0.4,\varphi = 0.5,\beta = 1.4,\eta = 0.2,\theta = 1.5,\nu = 1.8.$ The non-dominated solutions are obtained when using the multi-objective genetic algorithm are given in Figure 9, and the Pareto optimal solutions for various values are summarized in Table 2.

Figure 9. Pareto-front solution found by genetic algorithm.

Table 2. The Pareto optimal solutions for various values of λ.

From Figure 9 and Table 1, the three curves show that the increase in cost leads to a decrease in the waiting time of customers. ${F_p}$ increases as the average arrival rate increases. When ${T_p}$ approaches infinity, the limit value of the minimum expected cost ${F_p}$ can be regarded as the minimum cost.

Moreover, Figure 9 shows that the relationship between ${F_p}$ and ${T_p}$ is inversely related, approximating an exponential function with a negative exponent, and ${F_p}$ is positively related to λ. The original and adjusted coefficients of determination are ${R^2} = 99.6\% $ and ${R^2}(adj) = 99.6\% $, respectively, and the F-test of the regression equation has a p-value less than 0.01. Therefore, the regression model is appropriate. Figure 10 shows that the histogram of the regression residuals here is symmetric between the left and right sides, and the K-S normality test plot shows the scatter is essentially close to a straight line. The expected cumulative probabilities match well with the measured cumulative probabilities, indicating that the normal distribution is obeyed, the assumption of the regression is satisfied, and the regression relationship can be established. This regression model will help managers determine the minimum cost required to achieve a satisfactory level of service for a given estimate.

Figure 10. Histogram and K-S normality test for the regression residuals.

The least-square regression equation is established as follows:

(51)\begin{equation} {F_P} = - 8.039 + 8.846\lambda + 13.694T_P^{- 1} + 18.712{e^{- {T_p}}}. \end{equation}

6. Optimal strategy analysis

In this section, we give the reward-cost structure and focus on customers’ equilibrium joining strategies and the social benefit maximization problem. After service completion, each customer receives a reward R, and every customer pays the cost of remaining in the system, where the cost per unit of waiting in the system is C. All customers are risk-neutral and behave rationally to maximize their benefits.

Under the reward-cost structure as given above, the expected individual utility for the tagged customer who finds the server unavailable and decides to enter the orbit is

(52)\begin{equation} U(q) = R - CT(q), \end{equation}

where T(q) is given in (39).

When customers arrive at the system, they judge whether to join based on what they gain or lose. From the individual utility function, we can find the customers’ equilibrium joining probability ${q_e}$ as follows.

1. (i) If $U(1) \gt 0$, then the equilibrium strategy is ${q_e}=1$, this means that when customers choose to join the system, their expected individual utility is positive.

2. (ii) If $U(0) \lt 0$, then the equilibrium strategy is ${q_e}=0$, which means that when customers decide to join the system, their expected individual utility is negative.

3. (iii) In addition, the necessary and sufficient condition for ${q_e} \in (0,1)$ to be an equilibrium joining probability is that $U({q_e}) = 0$, this means that when customers choose to join the system, their expected individual utility is zero.

We continue with the problem of maximizing the social welfare per time unit.

The social welfare per time unit is

(53)\begin{equation} S(q) = {\lambda ^ * }R - CE(N), \end{equation}

where ${\lambda ^ * }$ is the effective arrival rate of customers, and

(54)\begin{align} \begin{split} {\lambda ^ * } &= \lambda {G_0}(1) + \lambda q[{G_1}(1) + {G_2}(1) + {G_3}(1) + {G_4}(1)]\\ &= \frac{{\lambda \mu {\kappa _3}{{\rm A}_1} + \lambda q(\beta + \varphi ){{\rm A}_2}}}{{{\kappa _2}[\theta \delta {\kappa _1} - \beta \mu \lambda {\kappa _3}]}}\pi (0,1). \end{split} \end{align}

From the above equation S(q), a socially optimal strategy ${q_{s}}$ is determined, which maximizes the social benefit per time unit.

Because the expressions obtained for individual and social benefits are very complex, it is hard to get specific results through traditional calculations. In the following analysis, we can use the particle swarm algorithm (PSO algorithm) to find the numerical solve this problem. When it comes to PSO algorithm, the most significant advantage is that it does not require too much analytic property of the objective function. It is an optimization algorithm based on swarm intelligence theory. During each iterative search, the particles in the swarm can dynamically adjust their position and velocity by tracking the two extremes of the swarm: the optimal solution P-best found by the particle itself and the optimal solution G-best found by the swarm. Through many iterations, the global optimal solution can be obtained.

In Figure 11, both ${q_e}$ and ${q_{s}}$ are monotonically decreasing concerning λ. Because more customers join the system, the orbital load will increase, which reduces the probability of customers joining the system. From Figure 12, both ${q_e}$ and ${q_{s}}$ decrease with respect to η. As η increases, customers’ waiting time on the orbit increases, which leads to a low probability of customers entering the system.

Figure 11. ${q_e}$ and ${q_{s}}$ versus λ $(\mu = 4.5,\beta = 0.3,\varphi = 0.1,\theta = 0.8,\eta = 0.1,\delta = 5,\nu = 10,R = 9,C = 1)$.

Figure 12. ${q_e}$ and ${q_{s}}$ versus η $(\lambda = 1.8,\mu = 5,\beta = 1.2,\varphi = 0.8,\theta = 1.5,\delta = 4,\nu = 5,R = 9,C = 1)$.

Figure 13 shows that the probabilities ${q_e}$ and ${q_{s}}$ increase with β. The repair rate increases, in this case, the server is quick to serve customers. In Figure 14, ${q_e}$ and ${q_{s}}$ increase with δ, delayed repair time becomes shorter, and the system performs repair services faster, so customers tend to join the system. As we expected, in Figure 15 ${q_e}$ and ${q_{s}}$ are monotonically decreasing concerning C. When the cost of waiting on the orbit increases, customers are unwilling to enter the system.

Figure 13. ${q_e}$ and ${q_{s}}$ versus β $(\lambda = 1.8,\mu = 4.5,\varphi = 0.8,\theta = 1.5,\delta = 2,\eta = 1,\nu = 7,R = 9,C = 1)$.

Figure 14. ${q_e}$ and ${q_{s}}$ versus δ $(\lambda = 1.8,\mu = 4.4,\beta = 2,\varphi = 0.5,\theta = 4.5,\eta = 1.3,\nu = 6,R = 9,C = 1)$.

Figure 15. ${q_e}$ and ${q_{s}}$ versus C $(\lambda = 2.5,\mu = 5,\beta = 0.3,,\varphi = 0.1,\theta = 0.7,,\eta = 0.1,\delta = 5,\nu = 10,R = 9)$.

Figure 16 shows that the maximum social benefit increases as the arrival rate λ and service rate µ increase. The reason is that the number of customers in the system increases, which contributes to increasing the social benefits. On the other hand, as the service rate increases, customers’ waiting times on the orbit decrease, so customers tend to join the system.

Figure 16. Maximum of social welfare versus λ for different values of µ $(\beta = 0.3,\varphi = 0.1,\theta = 0.8,\eta = 0.1,\delta = 5,\nu = 10,R = 9,C = 1)$ .

In Figures 17 and 18, the maximum social benefit increases as θ and β increase. Server availability increases due to the repair rate increased. However, the social benefit decreases as η and φ increase. As the breakdown rate increases, the server availability decreases. Thus, the social benefit decreases. It is shown in Figure 19 that the maximum social welfare increases as the reward R increases and decreases as the waiting cost C increases, which is consistent with our intuitive idea.

Figure 17. Maximum of social welfare versus θ for different values of η $(\lambda = 1.8,\mu = 5,\beta = 1.2,\varphi = 0.8,\delta = 4,\nu = 5,R = 9,C = 1)$ .

Figure 18. Maximum of social welfare versus β for different values of φ $(\lambda = 1.8,\mu = 4.5,\theta = 1.5,\eta = 1,\delta = 5,\nu = 7,R = 9,C = 1)$ .

Figure 19. Maximum of social welfare versus R for different values of C $(\lambda = 2.5,\mu = 5,\beta = 0.3,\varphi = 0.1,\theta = 0.7,\eta = 0.1,\delta = 4,\nu = 10)$ .

7. Conclusion

In this paper, we studied the M/M/1 retrial queueing system, which has passive breakdown with delayed repairs and active breakdown caused by negative customers during the busy period. Using the probability generating functions, we obtained the important performance measures, and we studied the effects of some parameters on the important performance measures of the model by numerical examples. We proposed a cost function to determine the optimal parameter settings for the system under stationary conditions. Moreover, we analyzed the equilibrium joining strategy and the socially optimal joining probability of customers and numerically analyzed the impact of some parameters on the maximum social welfare. A possible future research direction is to consider equilibrium strategies for unreliable M/M/C retrial queueing systems with negative customers, passive breakdown and delayed repairs as well as finding the number of servers such that the system incurs minimum costs.