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}

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.