key: cord-0061715-wnj1xcz8
authors: Yu, Miaomiao; Tang, Yinghui
title: Analysis of a renewal batch arrival queue with a fault-tolerant server using shift operator method
date: 2021-04-02
journal: Oper Res Int J
DOI: 10.1007/s12351-021-00635-4
sha: 3d8b3bb53941d0fd01f2432c174de647c617818d
doc_id: 61715
cord_uid: wnj1xcz8

Motivated by the fault tolerance for manufacturing, we investigate a renewal input bulk arrival queue with a fault-tolerant server, in which the server can keep working with a low service rate even if the partial failure occurs. Only when there are no customers in the system, the partial failure can be removed. To explore the performance measures of the queue, a more generic and simpler algorithm based on the right shift operator method for solving difference equations is employed to obtain the queue-length distributions at different time epochs. The significant feature of this algorithm lies in that it does not require the derivation of the transition probability matrix for the corresponding embedded Markov chain. Furthermore, we can resort to the queue-length distribution at the pre-arrival epoch to quickly get the expected sojourn time for an arbitrary customer. Finally, with the help of Pad[Formula: see text] approximation, several representative numerical examples are illustrated in tables and graphs, under which we show how to verify the correctness of our theoretical results through Little’s law.

A large number of papers in the field of queueing theory have appeared to analyze the waiting line system with an unreliable server. Among some earlier works in this area, Gaver (1962) studied the waiting line with interrupted service, including priorities. Avi-ltzhak and Naor (1963) and Thiruvengadam (1963) investigated some fundamental queueing problems with service facilities subject to breakdowns. Mitrany and Avi-Itzhak (1968) extended the analysis to the multiserver queue with service interruptions. Later, this model was revisited by Neuts and Lucantoni (1979) under the assumption that failed servers are repaired by one of the c repairmen. Over the past thirty years, authors like Sengupta (1990) , Takine and Sengupta (1997) , Tang (1997) , Li et al. (1997) , Madan (2003) , Ke (2003 Ke ( , 2006 , Wang (2004) , Gray et al. (2004) , Wang et al. (2010) , Choudhury and Deka (2012) , Jain et al. (2013) , Kumar et al. (2020) , Gorbunova and Lebedev (2020) , and Kumar et al. (2021) have studied some single arrival queueing systems with an unreliable server. A common assumption in the above literature is that as soon as the server fails, it instantaneously undergoes repairs. However, in real-life situations, due to some unavoidable reasons, the repair of the failed server may be delayed. This phenomenon stimulates Choudhury and Tadj (2009) , Choudhury and Kalita (2018) investigated the steady-state behavior of the unreliable M/G/1 queue with delay repair. On the other hand, from the mathematical and practical points of view, the case of batch arrival is more general, and also more challenging to handle. Many authors have contributed to the theory of batch arrival queues subject to unpredictable server breakdowns. Some notable works in this direction can be found in Ke and Lin (2006) , Huang (2010, 2012) , Choudhury and Tadj (2011), Singh et al. (2018) , Choudhury and Deka (2018) , Saggou et al. (2019) , Jain and Kaur (2020) . To get a more comprehensive summary of the recent research work on the topic of the unreliable queue, interested readers may refer to the review paper done by Jain et al. (2019a Jain et al. ( , 2019b .

It is worth noting that in the queueing systems mentioned above, the status of a server is usually modeled with two extreme states, working and failed. A working server is capable of serving customers at any instant of time. On the contrary, when the server fails, it is totally incapable of servicing and requires renewal or repair. However, such an assumption is not always true in reality. For example, a service facility or a manufacturing system built with fault tolerance capabilities will manage to keep operating (perhaps at a degraded level) in the presence of partial failure. Here, fault tolerance refers to the ability of a system to continue working without interruption when a partial failure occurs, wherein partial failures are often caused by the wear and tear of equipment components or the wrong action by an operator. The objective of creating a fault-tolerant system is to prevent disruptions arising in the service or manufacturing process, ensuring the high availability and business continuity of a stochastic service system. To theoretically analyze such phenomenon, Kalidass and Kasturi (2012) developed a new concept of working breakdowns in queueing theory. In a working breakdown queue, the server works at a lower service rate rather than completely stopping service during the breakdown period. Many researchers motivated by their work have extended this type of queue in different frameworks. Kim and Lee (2014) further studied an M/G/1 queue with disasters and working breakdowns, in which the system is equipped with a substitute server for providing the working breakdown services to arriving customers. Liu and Song (2014) extended the idea in Kalidass and Kasturi (2012) to the Markovian batch arrival queue. Applying the matrix-geometric method (Neuts 1981) , Liou (2015) also investigated M/M/1 unreliable queue subject to working breakdowns and impatient customers. Other recent works on queues with working breakdowns can be found in Li et al. (2013) , Yen et al. (2016) , Chen (2018), Ye and Liu (2018) , Jain et al. (2019a Jain et al. ( , 2019b , Xin (2019), and Gao et al. (2019) . It is, however, slightly regrettable that very few authors have studied the working breakdown queue under general renewal arrivals. To date, only a small amount of literature has covered this topic. Using the embedded Markov chain and matrix analytic approach, Jiang and Liu (2017) considered the GI/M/1 queue with disasters and working breakdowns in a multi-phase service environment. Utilizing the supplementary variable technique (Cox 1955) , Yang and Cho (2019) presented a recursive algorithm for computing the stationary queue-length distribution in N-policy GI/M/1/K queue with modified working breakdowns. Here, Yang and Cho did not assume that service and maintenance can be performed simultaneously. In their model assumptions, the repairman repairs the partially failed service facility when there are no customers queueing up for service. We think this assumption may be valid in some cases. For example, as the COVID-19 pandemic has spread across the globe, the severe shortage of medical masks caused by the health crisis quickly becomes a significant issue of the pandemic. To address the growing demand for masks during the outbreak of coronavirus disease, Honeywell, a leading manufacturer of personal protective equipment, quickly ramps up production by reducing regular equipment maintenance frequency. Such maintenance schedule adjustment can be regarded as a fault-tolerance mechanism in the production process that simultaneously provides a timely response to customer needs. This example also gives us a realistic background to study the queueing system with fault-tolerant operations.

It is an indisputable fact that queueing models with customers arriving in batches rather than singly have many applications in practice, for example, in flexible manufacturing systems. However, based on the above literature review, we may see that current research has not paid much attention to the renewal batch arrival queue with working breakdowns or with fault tolerance characteristics. Such a situation prompted us to study a bulk queue with renewal input and fault-tolerant operating modes. Also, to simplify the analytical study of the model, an alternative yet simple method with resorting to the supplementary variable technique and the shift operator technique in solving difference equations is employed to analyze our model. We note that in the past one year, Barbhuiya and Gupta (2019a , 2019b used this unique method to reconsider the algorithm for computing the queue-length distributions in some renewal input bulk arrival queues. From their pioneering work, it was revealed that the algorithm based on the combination of these two techniques is easy to understand, and can be conveniently implemented by applying a suitable software package such as Maple, Mathematica or Matlab. More importantly, this method requires neither derivation of the transition probability matrix of the Markov chain embedded at arrival instants nor the inversion of any probability generating function. We quickly realized that a more generic version of this method could be developed and applied to analyze our queueing model presented in this work.

The remaining part of our paper is organized as follows. We first describe the system and provide basic assumptions in Sect. 2. The differential-difference equations governing the queueing model are framed by the supplementary variable technique in Sect. 3. In Sect. 4, with the help of the theory of difference equation, we provide the procedure to analyze the steady-state queue-length distributions at different epochs of our model. Moreover, using the queue-length distribution at the pre-arrival epoch, we also evaluate the sojourn time of an arbitrary customer in Sect. 5. To validate the correctness of our theoretical results, we present some numerical examples in Sect. 6. Section 7 concludes the paper and directs possible further studies.

The mathematical model for this work describes as follows:

• Consider GI X /M/1 queue with a fault-tolerant server wherein batches of customers arrive at epochs 0 = 0 , 1 , 2 , … , n , … . The number of arrivals at each epoch is given by a random variable X having general distribution g i = Pr {X = i} , i = 1, 2, … , and probability generating function

In most real-life situations, batch sizes are not infinite but finite. Thus, from a realistic and computational point of view, we suppose that the maximum batch size is b throughout this paper. As a result of this assumption, the p.g.f. and mean of the random variable X are given by

• The inter-arrival times n+1 − n > 0 , n = 0, 1, 2, … , are independent and identically distributed (i.i.d.) random variables with common distribution function A(t) and probability density function (p.d.f.) a(t). Let the Laplace-Stieltjes transform of this distribution be denoted by a * (s) = ∫ ∞ 0 e −st dA(t) and let the mean interarrival time be denoted by 1∕ , where 0 < 1∕ = − d ds a * (s) | | |s=0 < ∞. • The customers are served individually by a single server on a first-come, firstserved basis. If S n is the service time of the nth customer to be served in normal state, then it is assumed that S n , n = 1, 2, 3, … is a sequence of positive i.i.d. random variables with the common exponential p.d.f. 0 e − 0 t , t > 0. • A partial failure of the server could occur during the normal busy period, and the time until the random partial breakdown of the server, denoted by L, is assumed to be exponentially distributed with rate . • When a partial failure occurs, the server continues working in degraded mode at a reduced service rate rather than entirely halting service. Let S n represent the service time of the nth customer to be served in defective state, and we also assume that the sequence S n , n = 1, 2, 3, … are i.i.d. positive random variables following exponential distribution with parameter 1 ( 1 < 0 ).

• Given customers' sensitivity and emotional reaction to delay, and also to ensure continuity and efficient service process, the defective server is repaired after the system becomes empty. In other words, when there are customers present in the queue, the server cannot be repaired, even if the server is undergoing repairs. The repair time of the defective server is exponentially distributed with a mean rate . • It is further supposed that various stochastic processes involved in the queueing system are mutually independent. Additionally, for the stationary analysis of the model, we demands that =ḡ ∕ 1 < 1 (see the proof of Theorems 1 and 2).

Let N(t) and (t) indicate the number of customers in the system (including the one being served) and the state of the server at time t, respectively. Here, (t) is a binary random variable. (t) = 0 represents that the server is in a normal state, while (t) = 1 represents that the server is defective. Moreover, the supplementary variable R(t) corresponding to the remaining inter-batch arrival time at time t is used, which provides Chapman-Kolmogorov equations governing our current model. Thus, the state of the system at time t could be described by a multivariate stochastic process

In this process, the non-Markovian queueing process becomes Markovian by having the necessary information so that the future of the process depends only on its current state. To ease us into the analysis of such a queue, we define the joint probabilities as follows

In steady-state, let us further define

To obtain the steady-state probabilities of the queue size Q n,0 and Q n,1 , we first construct the differential-difference equations for describing the evolution of the model by observing the state of the queue at two consecutive time epochs t and t + Δt . Using the probabilistic argument, we can derive the following set of equations: (

Since the Laplace transform are extremely useful in the solution of complicated differential-difference equations presented above, we define the Laplace transform of Q n,m (x) as Meanwhile, we also notice that Q * n,m (0) = Q n,m . Then, taking the Laplace transform on both sides of Eqs. (1)-(6) and using the relation − ∫ ∞ 0 e −sx d dx Q n,m (x)dx = Q n,m (0) − sQ * n,m (s) , we thus obtain the transformed equations as below

Adding Eqs. (7)-(12), it yields

Taking the limit as s → 0 , and using L'Hôspital's rule one time in the above Eq. (13), we have

Here, the normalization condition implies that Substituting Eq. (15) into Eq. (14), we get Employing Eq. (16), we may show that if Q − n,m represents the steady-state probability of having n (n ≥ 0) customers in the system and the server being in state m (m = 0, 1) immediately before batch arrival, then Q − n,m = 1 Q n,m (0) . To show this, it should be noted that Q − n,m is proportional to Q n,m (0) and

Now, we have already seen that once Q n,m (0) and Q * n,m (s) are obtained, we can further use them to get the queue-length distributions at the instant of a batch arrival

( Q − n,m ) and arbitrary epoch ( Q n,m ). We will address how to deal with this problem in the next section.

For analysis purposes, the discrete variable n will be considered to be the independent variable, while the function values Q n,m (0) and Q * n,m (s) (m = 0, 1) will be the dependent variables. For the sequences Q n,m (0), n ≥ 0 and Q * n,m (s), n ≥ 0 , we define a right shift operator on the above two sequences and set j Q n,m (0) = Q n+j,m (0) and j Q * n,m (s) = Q * n+j,m (s) , j ≥ 1 . Thus, Eq. (9) can now be written in terms of operator as follows By setting n ≥ 0 instead of n ≥ b and substituting s = + 0 − 0 , Eq. (18) reduces to the below homogeneous difference equation with constant coefficients According to the fundamental theory of difference equation, the corresponding auxiliary equation is said to be the characteristic equation of Eq. (19). Next, we will apply Rouché 's theorem to find the number of roots of Eq. (20).

is holomorphic inside and on the closed contour |z| = 1 + . Thus, according to the Taylor's theorem for analytic complex function there exists a power series

and ℂ is any closed contour around 1 and lying completely inside |z| ≤ 1 + . Then, employing the Taylor expansion for H(z), we may estimate |f(z)| and |g(z)| on the simple closed curve |z| = 1 − , where > 0 and is also a sufficiently small quantity. Thus, on

Consequently, as all the conditions of Rouché 's theorem (see Klimenok 2001) are satisfied, f (z) + g(z) has exactly b zeros inside the unit circle, since f(z) has b. ◻ It has been shown by many authors that in queueing theory, roots of the characteristic equation are well structured, and they are generally distinct, see, e.g., Tijms (2003) and Chaudhry et al. (1990) . Therefore, if we assume that the roots of Eq. (20) are distinct, denoted by r j , j = 1, 2, … , b , the general solution of Eq. (19) can be expressed in the form where j (j = 1, … , b) are real or complex constants to be determined. Substituting Eq. (21) into Eq. (18) yields It is now necessary to turn to determine the unknown function Q * n,0 (s) from the above equation. Treating s as a fixed constant, Eq. (22) 

Thus, a particular solution of the nonhomogeneous Eq. (22) can be given as follows

Adding this particular solution to the general solution of the associated homogeneous equation, the solution of Eq. (22) is presented as

We now find the conditions under which Q * n,0 (s) has the same expression as in Eq. (25) for 1 ≤ n ≤ b − 1 . To this end, substituting Eq. (21) into Eqs. (8) and (9), and comparing the last term of the right hand side of the above two equations, we may observe that the unknown constants j satisfy the relationship Respectively setting n = b − 1, b − 2, … , 1 in Eq. (26) and noting that g b ≠ 0 , we can construct a system of b − 1 linear equations in b variables which can be employed to determine the unknown constants j (j = 1, 2, … , b) in later analysis

If the above results hold true, for any n ≥ 1

Next, let us once again apply the right shift operator to Eq. (12). Then, we have By setting s = 1 − 1 in Eq. (29) and substituting Eq. (28) into Eq. (29) yields Adopting Rouché 's theorem similar to that used for the previous Theorem 1, the following Theorem 2 shows that under certain conditions, the characteristic equation of the above difference equation also has exactly b roots inside the unit circle |z| = 1 . We shall use the letters 1 , 2 , … , b to denote these b roots.

has exactly b roots inside the unit disk.

Furthermore, combining the results of Theorems 1 and Theorem 2, the obvious stability condition of the current queueing system is that ḡ 1 < 1 . By analogy with the corresponding procedure to compute Q * n,0 (s) , the general solution of the nonhomogeneous Eq. (30) is given by

Here, the first term in the right-hand side of Eq. (31) is a solution to the associated homogeneous equation of Eq. (30), and k 1 , k 2 , … , k b are the arbitrary constants that can be obtained from later analysis. On the other hand, the second term is a particular solution of Eq. (30). Substituting Eqs. (28) and (31) into the right-hand side of Eq. (29), we have

which is also a nonhomogeneous difference equation with constant coefficients. Similar to the discussions surrounding Eq. (22) (s) . Summing over all permissible n from b to ∞ and taking the limit as s → 0 , the formula ∑ ∞ n=b Q * n,1 (0) = ∑ ∞ n=b Q n,1 ≤ 1 clearly holds. This means that the undetermined constant D 2 = 0 . Otherwise, the limit of

as s approaches 0, is infinity. Therefore, the solution of Eq. (32) takes the below form (32)

We now find the condition under which the expression for Q * n,1 (s) presented in Eq. (34) also holds when 1 ≤ n ≤ b − 1 . Insert Eqs. (28) and (32) into Eqs. (11) and (12), respectively. By comparing the last term of the right hand side of Eqs. (11) and (12), we need to imposes the following restrictions on the constants k j and j .

By putting n = b − 1, b − 2, … , 1 in Eq. (35) and noting that g b ≠ 0 , Eq. (35) can be conveniently written in the form of a linear equation with 2b variables, i.e., k 1 , k 2 , … ,

In other words, for 1 ≤ n ≤ b − 1 and n ≥ b , if the above Eq. (36) holds, Q * n,1 (s) has the uniform expression as follows:

Setting s = and using Eqs. (31) and (37) in Eq. (10), we can obtain another linear equation in the same 2b variables (34)

By summing over all possible n from 0 to ∞ in Eq. (16) and employing Eqs. (21) and (31), we can further derive a relationship for j and k j , j = 1, 2, … , b.

So far, by the use of Eqs. (27), (36), (38) and (39), we have established a system of 2b linear equations in 2b variables which can be solved to obtain the constants j 's and k j 's, j = 1, 2, … , b . Once the arbitrary constants are determined by solving the linear algebraic equations, the stationary queue-length distributions at pre-arrival and arbitrary epochs are given, respectively, by

Remark 1 Our queueing model is a generalization of the classical M/M/1 queue. More precisely, such a model reduces to M/M/1 queue when → 0 , → ∞ , 0 = 1 , g 1 = 1 and a * (s) = s+ . For this case b = 1 , thus the single characteristic root inside the unit circle is r 1 = 1 = ∕ 1 . Hence, we may compute the probability that an arriving customer finds an empty system (denoted as Q − 0 ) from Eq. (40) and obtain

On the other hand, substituting r 1 = 1 = ∕ 1 into Eq. (39), we can find that 1 +k 1 satisfies the relationship 1 +k 1 = 1 − 1 = (1 − ) = Q − 0 . The relationship also implies that Q − n = Q − n,0 + Q − n,1 = 1 +k 1 r n 1 = (1 − ) n , n = 1, 2, … , where Q − n represents the probability that there are n customers in the system just before the arrival of a customer. Therefore, according to the PASTA property (Poisson Arrivals See Time Averages), we get exactly the same results that have been reported in the existing literature on queueing theory (see Gross and Harris 1985) .

(40)

Furthermore, setting s = 0 in Eq. (10), we obtain

With the help of the normalizing condition ∑ ∞ n=0 Q n,0 + ∑ ∞ n=0 Q n,1 = 1 , we may further obtain This completes the analysis of queue-length distribution at different epochs.

The sojourn time T A is the total time that an arbitrary customer in an arriving batch spends in the system until it departs from the system. Using the pre-arrival epoch probabilities Q − n,0 and Q − n,1 that we have derived in Sect. 4, we will investigate the Laplace-Stieltjes transform and the expectation of an arbitrary customer's sojourn time in this section. For this purpose, we introduce a tagged customer who may arrive at a random position in the queue. Clearly, it may be seen that a tagged customer arrival may belong to one of the following cases:

A batch containing the tagged customer arrives at the system during the normal working period and finds n customers already present in the system. Meanwhile, the number of customers that arrive in the same batch as the tagged customer, but enter service before the tagged customer is l (l = 0, 1, … , b − 1) . To discuss the sojourn time of the tagged customer, this case can further be divided into two subcases: (i) The time until the random partial breakdown of the server is no less than the total service time of n + l + 1 customers. That is to say, the tagged customer completes its service and leaves the system in the normal working period; (ii) There are only m (m = 0, 1, … , n + l) service completions during the normal working period. In other words, the tagged customer leaves the system when the server undergoes in the defective state.

Case 2. An arriving batch sees the state of the system is (n, 1), n ≥ 0 , and the (l + 1) th customer in this batch is the tagged customer. Then, the sojourn time of the tagged customer equals the total customers' service times ahead of him plus his service time. Since the server is defective, these services are only provided with the low service rate 1 .

Let W A (t) be defined as the probability distribution function of T A , and g − l denotes the probability of l number of customers ahead of a randomly selected tagged customer within the batch. From Burke (1975) and Chaudhry and Templeton (1983) , g − l is given by g − l = 1 g ∑ ∞ j=l+1 g j . Based on the above different cases, we have, from the theorem of total probability, The Laplace-Stieltjes transform of the distribution of T A is given by

Hence, the first and second moments of T A may be found from the Laplace-Stieltjes transform as

From the numerical examples presented in the next section, we will see that Eq. (46) can provide us an effective way to validate the correctness of our theoretical analysis results. Moreover, using Eqs. (46) and (47) we can compute the variance of T A as

To demonstrate the working schemes of the difference equation approach based on the right shift operator, we first describe the solution algorithm for calculating the steady-state probabilities Q − n,0 , Q − n,1 , Q n,0 and Q n,1 , for n ≥ 0 . Given the values of 0 , 1 , , , the Laplace-Stieltjes transform expression of the inter-batch arrival time a * (s) and the probability mass function g i = Pr {X = i} , i = 1, 2, … , b , the steps of the solution algorithm are stated as follows:

• Step 1: Find the roots of the following two characteristic equations (42) and (43), compute Q n,0 and Q n,1 , for n ≥ 0.

All the calculations are performed on a PC having Corei7 processor at 3.20 gigahertz with 16 gigabytes RAM using Mathematica and Matlab software packages. We use Mathematica software to find the roots of the associated characteristic equations, and then write a Matlab code to solve a linear system of equations in 2b variables. Next, to illustrate the solution algorithm, we provide three numerical examples where the inter-batch arrival time distributions are 2-stage Erlang, uniform and deterministic, respectively. A variety of numerical results have been presented in self-explanatory tables and graphs. The notations used in these tables are the same as those defined earlier in the previous sections except E T A Little , which denotes the average sojourn time of an arbitrary customer evaluated through Little's formula.

Example 1 Consider the E X 2 /M/1 queueing system with fault tolerance capabilities. The 2-stage Erlang distribution is made up of two independent and identical exponential stages, each with mean 2.5. In this case, we have a * (s) = 0.4 ( i = 1, 2, … , 10 ) with mean value ḡ = 1.814776 . For computation purpose, we fix the other parameters as 0 = 0.7 , Table 3 , an effective approach is provided to verify the correctness of our numerical as well as analytical results. Actually, there are two different ways to get the average sojourn time for an arbitrary customer. One way is to substitute Q − n,0 and Q − n,1 into Eq. (46) and calculate E T A directly. The other is to apply Little's formula to obtain E T A = L s ∕ḡ , where L s = ∑ ∞ n=1 n(Q n,0 + Q n,1 ) . The results from numerical computation (see the bottom of Table 3) indicate that E T A evaluated through Eq. (46) exactly matches with the one obtained from Little's formula. It also implies that the theoretical analysis and numerical experiments performed in this paper is valid and reliable.

In the first example, we consider the case when inter-arrival times of groups have Erlangian distribution of order 2. We note that such an arrival process is a particular case of the batch Markovian arrival process (BMAP). The corresponding numerical results can be much easier obtained by using the well-known matrix geometric method. To further demonstrate the universality of the proposed method, some other examples with an inter-batch arrival time that does not belong to the class of PH distribution are given below.

Example 2 Consider the U X /M/1 queue with fault tolerance capabilities, where the inter-batch arrival times are independent random variables, each distributed uniformly on interval 0 to 2∕ . If we set = 0.2 , the uniform distribution has the Laplace-Stieltjes transform a * (s) = 1 10s (1 − e −10s ) . The maximum batch size is b = 12 , and take the probability mass function of batch size to be

e −2 2 n n! , i = 1, 2, … , 12 with mean value ḡ = 1.930825 . Since a * (s) is a transcendental function, the corresponding characteristic equations cannot be directly solved by using the standard Mathematica commands. For the purpose of finding the roots of the characteristic equations, we wish to approximate a * (s) by means of a rational approximation ℙ 1 (s)∕ℙ 2 (s) , where ℙ 1 (s) and ℙ 2 (s) are polynomials of degree m 1 and m 2 respectively. It is well known that the Padé approximation is a particular and classical type of rational fraction approximation. Many practical applications have proven that it is the best approximation of a function by a rational function of a given order. In this example, we shall use the so-called Padé rational approximation of degree (15, 16) to approximate a * (s)

Here, calculations with Padé approximant of a * (s) are straightforward and can be performed with Mathematica command "PadeApproximant". How the choice of Padé ( m 1 , m 2 ) affects the accuracy of numerical results and how one can trade off between computation-time and accuracy is discussed in detail in the work of Singh (2014) . By the fixed values of certain parameters 0 , 1 , and as 0.75, 0.5, 0.003 and 0.1 respectively, we present the numerical results for the roots of the characteristic equations inside the unit circle in Table 4 .

Substituting these roots into Eqs. (27), (36), (38) and (39), and solving a simultaneous set of twenty-four linear equations, Table 5 gives the numercal results of j and k j for j = 1, 2, … , 12.

The summary of the calculations for the queue-length distributions at prearrival and arbitrary epochs is shown in Table 6 . At the same time, the expected sojourn time for an arbitrary customer estimated in two different ways is also summarized at the bottom of Table 6 . We may see that the numerical experiments verify our theoretical results and show their correctness once again.

Example 3 Consider the D X /M/1 queue with a fault-tolerant server, in which the batch arrivals are equally spaced in time. For this inter-batch arrival time distribution a * (s) = e − s , here we take = 0.2 . Now, let us conduct the numerical experiment with the following parameters: 0 = 0.8 , 1 = 0.6 , = 0.0025 and = 0.15 . Additionally, the corresponding numerical results were obtained by assuming that the batch size distribution has a p.g.f. G(z) = 0.25z + 0.5z 2 + 0.1z 4 + 0.1z 6 + 0.05z 8 . This suggests that the batch size can be either 1, 2, 4, 6, or 8 with 25% , 50% , 10% , 10% and 5% probabilities, respectively. Since the deterministic inter-batch arrival time does not have a rational Laplace-Stieltjes transform, Mathematica software package cannot solve the characteristic equations (transcendental equations) directly. Through the Padé approximation of degree (8, 9) , we also approximate the Laplace-Stieltjes transform of the inter-batch arrival time distribution e −5s with a rational function of the type ℙ 1 (s)∕ℙ 2 (s):

Letting s = 0.8025 − 0.8z and s = 0.6 − 0.6z in the above expression, respectively, and putting them into the characteristic equations, Table 7 gives all distinct roots that are found within the contour of a unit circle |z| = 1. The unknown constants j and k j for j = 1, 2, … , 8 can be found in the same manner as described earlier. (42) and (43), we compute the probability distributions of the queue length at two different epochs, and the numerical results are listed in Table 9 . It is also found that E T A matches exactly with the mean sojourn time calculated using Little's law.

In this example, the computer software, e.g., Mathematica and Matlab, are used to compare three configurations in terms of their E T A and L s for three different inter-batch arrival time distributions: 2-stage Erlang, uniform and deterministic. We first perform a comparison for the average sojourn time of an arbitrary Table 2 The values of j and k j for j = 1, 2, … , 10 in Example 1 customer using the assumption that these distributions have the same mean 5 but different standard deviations. We choose 0 = 0.9 , 1 = 0.75 , = 0.0025 , and suppose X is a discrete uniform random variable on the consecutive integers 1, 2, … , 6, so that the traffic intensity = 0.883333 . At the same time, to evaluate the impact of the repair rate on the value of E T A , we vary the values of from 0.05 to 0.25 and draw the plot of E T A as a function of . Figure 1(a) depicts that E T A decreases as increases. That is to say, moderate shortening the sojourn time of an arbitrary customer can be achieved by choosing a higher repair rate. On the other hand, we also observe that the expected sojourn time is larger for 2-stage Erlang distribution as compared to uniform and deterministic distributions. We think in most of the cases, by comparing E T A in terms of three different distributions of the inter-batch arrival time usually yields E T A E 2 > E T A U > E T A D . This is due to the fact that the standard deviation of the deterministic distribution equals zero. In Figure 1(b) , under the same parameter settings, it can be seen that the average queue length in the case of 2-stage Erlang inter-batch arrival time distribution with a higher standard deviation is larger than the one with lower standard deviations. Again deterministic distribution with zero standard deviation yields the lowest average queue length. Such results indicate that the inter-batch arrival time plays a significant role in determining the performance of the queueing system.

In this example, we perform a comparative analysis on E[T A ] and L s based on changes in specific random arrival batch size distributions and assume that the inter-batch arrival times have a common 3-phase hyperexponential distribution, whose density and Laplace transform are given by A k-phase hyperexponential distribution is frequently used in queueing theory to model the distribution of the superposition of k independent events. For instance, the arrival of different types of customer to a single queueing station is often modeled as a hyperexponential distribution. To demonstrate the performance evaluation, we consider the following three batch size distributions:

Case 1: Batch size X follows the discrete uniform distribution on the interval [1, 2, … , 16].

Case 2: The batch size distribution has the probability mass function: g 4 = 0.25 , g 8 = 0.25 , g 10 = 0.25 , and g 12 = 0.25.

Case 3: The arriving batch size follows a 6-8-12 distribution with probability mass function g 6 = 0.25 , g 8 = 0.5 , and g 12 = 0.25.

All these distributions are characterized by the average batch size equal to 8.5. However, these are qualitatively different in that they have different variances, the variances of above probability distributions are 21.25, 8.75 and 4.75, respectively. The numerical experiment has been done by fixing other system conclude that the variance of batch size has obvious impact on performance metrics like average system size and average sojourn time.

In order to mitigate the severe potential consequences of system failure, fault-tolerant operating has become a critical attribute of the manufacturing system. In this paper, we have developed a GI X /M/1 queue with a fault-tolerant server to study the behavior of key performance indicators for such a system. The mathematical analysis of the queueing model is carried out by employing a combination of two methods. One is the well-known supplementary variable technique, and the other is the shift operator method in difference equations theory. The most important reason why we adopt the above approach is that we want to avoid discussing the transition probability matrix associated with the embedded Markov chain. Because with the increase in the complexity involved in the queue, the derivation of the transition probability matrix becomes extremely difficult, and many practitioners cannot fully master these complicated and tedious probabilistic arguments. In essence, it is our intention to solve a practical queueing problem in the most simple and pragmatic way using the roots of the model's characteristic equation. From the algorithmic steps and numerical experiments presented in our work, we may also see that the methodology is nicely implemented in flexible scientific computing software that is easy to use even for users who are not professional mathematicians. A feasible plan for further study is to consider the discrete-time counterparts of this queue. Since several possible events may occur simultaneously in discrete-time queues, the analysis and computation will become more elaborate under such a situation. 

Some queueing problems with the service station subject to breakdowns

A difference equation approach for analysing a batch service queue with the batch renewal arrival process

Discrete-time queue with batch renewal input and random serving capacity rule: GI X /Geo Y /1. Queueing Syst

Analytical and computational aspects of the infinite buffer single server N policy queue with batch renewal input

Delays in single-server queues with batch input

Robustness of rootfinding in single-server queueing models

System reliability analysis of retrial machine repair systems with warm standbys and a single server of working breakdown and recovery policy

A single server queueing system with two phases of service subject to server breakdown and Bernoulli vacation

A batch arrival unreliable server delaying repair queue with two phases of service and Bernoulli vacation under multiple vacation policy

An M/G/1 queue with two types of general heterogeneous service and optional repeated service subject to server's breakdown and delayed repair

An M/G/1 queue with two phases of service subject to the server breakdown and delayed repair

The optimal control of an M X /G/1 unreliable server queue with two phases of service and Bernoulli vacation schedule

The analysis of non-Markovian stochastic processes by the inclusion of supplementary variables

Analysis of a discrete-time repairable queue with disasters and working breakdowns

A waiting line with interrupted service including priorities

Queueing system with two input flows, preemptive priority, and stochastic dropping

A queueing model with multiple types of server breakdowns

Supplementary variable technique (SVT) for non-Markovian single server queue with service interruption (QSI)

Policy for unreliable server bulk queue with Bernoulli feedback

Unreliable server M/G/1 queue with multi-optional services and multi-optional vacations

Performance modeling of fault-tolerant machining system with working vacation and working breakdown. Arab

The GI/M/1 queue in a multi-phase service environment with disasters and working breakdowns

Computational analysis of the queue with working breakdowns and delaying repair under a Bernoulli-schedule-controlled policy

A queue with working breakdowns

queueing system with server vacations, startup and breakdowns

queue under hysteretic vacation policy with an early startup and un-reliable server

Analysis of an unreliable server

Analysis of batch arrival queue with randomized vacation policy and an unreliable server

Maximum entropy solutions for batch arrival queue with an un-reliable server and delaying vacations

The M/G/1 queue with disasters and working breakdowns

On the modification of Rouche's theorem for the queueing theory problems. Queueing Syst

Performance analysis of multi-processor twostage tandem call center retrial queues with non-reliable processors

Performance analysis of an unreliable M/G/1 retrial queue with two-way communication

Reliability analysis of M/G/1 queueing systems with server breakdowns and vacations

Equilibrium customer strategies in Markovian queues with partial breakdowns

Markovian queue optimisation analysis with an unreliable server subject to working breakdowns and impatient customers

The M X /M/1 queue with working breakdown

An M/G/1 type queue with time-homogeneous breakdowns and deterministic repair times

A many-server queue with service interruptions

A Markovian queue with N servers subject to breakdowns and repairs

The analysis of unreliable M [X] /G/1 queuing system with loss, vacation and two delays of verification

A queue with service interruptions in an alternating random environment

Performance analysis of bulk arrival queue with balking, optional service, delayed repair and multi-phase repair

Analysis of queueing time distributions for MAP/D N /1 queue

A single server queue with service interruptions

A single server M/G/1 queueing system subject to breakdowns-some reliability and queueing problems

Queuing with breakdowns

A First Course in Stochastic Models

queue with second optional service and server breakdowns

Comparison of two randomized policy M/G/1 queues with second optional service, server breakdown and startup

Analysis of the N-Policy GI/M/1/K queueing systems with working breakdowns and repair

Analysis of MAP/M/1 queue with working breakdowns

Reliability and sensitivity analysis of the controllable repair system with warm standbys and working breakdown

The authors would like to thank the anonymous reviewers for their valuable comments and suggestions to improve the quality of this paper. This research was supported by the National Natural Science Foundation of China (Grant No. 71571127), and the National Office for Philosophy and Social Sciences of P. R. China (Grant No. 20BGL109).

 Table 4 The roots of the characteristic equations with modulus less than one in Example 2 j r j j r j j j j j 1 0.706120 7 0.157883 + 0.199331 1 0.939953 7 −0.144906 + 0.215858