Steady-State Analysis of a Flexible Markovian Queue with Server Breakdowns

A flexible single-server queueing system is considered in this paper. The server adapts to the system size by using a strategy where the service provided can be either single or bulk depending on some threshold level c. If the number of customers in the system is less than c, then the server provides service to one customer at a time. If the number of customers in the system is greater than or equal to c, then the server provides service to a group of c customers. The service times are exponential and the service rates of single and bulk service are different. While providing service to either a single or a group of customers, the server may break down and goes through a repair phase. The breakdowns follow a Poisson distribution and the breakdown rates during single and bulk service are different. Also, repair times are exponential and repair rates during single and bulk service are different. The probability generating function and linear operator approaches are used to derive the system size steady-state probabilities.


Introduction
Markovian queueing models have Poisson arrivals and exponential service times. Because they are (arguably) easy to analyze, they are often used as a first step in analyzing more difficult queueing systems. They also yield practical results that present no difficulty in implementation. Early models were Markovian and found application in the telephone industry. However, since then, they also found other areas of application such as computer, transportation, and production systems.
The literature on Markovian queueing systems is huge. To cite a few of the latest, Wang et al. [1] use a Markovian queue to model a passenger-taxi system. They design a social benefit function and look for the system parameters that optimize the system operations.
Jain et al. [2] investigate a theoretical model of a Markovian system where the server takes vacations. During a vacation, the server does not stop serving customers but reduces his service rate. Also, customers can get discouraged and may not join the queue. Jain et al. study the transient behavior of the system using probability generating functions.
While researchers use Shannon entropy to measure randomness in queueing systems, Srivastava [3] uses Renyi's measure of entropy to quantify uncertainty in Markovian queueing systems with finite and infinite capacity.
Estimation of the parameters of Markovian queueing systems using statistical techniques and simulation is also reported in numerous papers, see for example Refs. [4][5][6].

Model Formulation
Oftentimes, a service or manufacturing firm can process its customers either singly or in batches, see case study below. It may not be economical to process customers singly, and it may be impossible to always process them in batches. In this case, a service discipline where the firm decides to process customers singly at times and in batches at other times is more appropriate. We therefore consider in this paper a single-server Markovian queueing system where the switch from one service discipline to the other is triggered by some constant integer c ≥ 2. If the number of customers in this system is less than c, then the server processes customers one at a time (single mode). If the number of customers in greater than or equal to c, then the server provides service simultaneously to a group of c of customers (batch mode).
To define the rest of the notation, we assume that customer arrivals follow a Poisson process with positive rate λ. In single mode (i.e., when the number of customers in the system is less than c), service follows an exponential distribution with positive rate µ 1 . In batch mode (i.e., when the number of customers in the system is greater than c), service follows an exponential distribution with positive rate µ 2 > µ 1 . To mark the transition from single mode to batch mode, we assume that when the number of customers in the system is equal to c, then service is batch but with rate µ 2 − µ 1 . We also assume that the server is unreliable and may break down, either in single mode or in batch mode. Breakdowns occur according to a Poisson process with positive rate α 1 in single mode and α 2 in batch mode. A breakdown is followed by a repair of the server and repair times follow an exponential distribution with positive rate β 1 in single mode and β 2 in batch mode. When the server is unavailable, customers are allowed to join the queue in single mode, but not in batch mode, to avoid large queue lengths.
Let X(t) denote the number of customers in the system at time t and P n (t), n = 0, 1, 2, · · · denote the probability of n customers in the system at time t. The process {X(t); t ≥ 0} is a continuous-time Markov chain, and the corresponding rate-transition diagram is depicted in Figure 1.
Since the server may be either working or down, we introduce W n (t), n = 0, 1, 2, · · · the probability of n customers in the system at time t when the server is in a working state, and F n (t), n = 0, 1, 2, · · · the probability of n customers in the system at time t when the server is in a failing state. Note that we readily have P n (t) = W n (t) + F n (t), the probability of n customers in the system at time t, regardless of the server state. Writing the Chapman-Kolmogorov equations in the case where the server is in a working state, Similarly, writing the Chapman-Kolmogorov equations in the case where the server is in a failing state, we have Taking the limit in both sets of difference-differential equations as t → ∞ yields the balance equations By substitution of (10), (11), and (12) into (7), (8), and (9), respectively, we obtain (λ + µ 2 )W n = λW n−1 + µ 2 W n+c , n ≥ c.
These difference equations are solved in the next section.

Model Solution
We use two methods to solve the set of Equations (13)- (15). The first method involves probability generating functions (PGF). The second one involves the concept of linear operators. Our aim is to solve the considered problem using the above two different methods and then to compare the obtained probabilities.

Analytical Method Using PGF
The procedure is to find a closed-form expression for the PGF. Then, if possible, expand it as a power series. If not, the probabilities are obtained through successive differentiations.

Computations of the Probabilities P n
In this method, to find the steady state probabilities W n , F n , and P n , we introduce for |z| ≤ 1 the probability generating functions: To start, we define We multiply, for 1 ≤ n ≤ c − 1, both sides of Equation (14) by z n+c : Taking the summation of these equations over n from 1 to c − 1 yields: which can be rewritten Now, we multiply, for n ≥ c, both sides of Equation (15) by z n+c : Similarly to the previous case, we take the summation of these equations over n from c to ∞ to get the following: which is equivalent to Rearranging terms, we have From Equation (19) we have Now replace µ 2 S 2 (z) with its expression in (23) to obtain: Simplifying the RHS of this expression then solving for W(z), we get the PGF of the working steady-state probabilities Finally, using (7) we get the following where (11) and (12), we obtain the PGF of the failing steady-state probabilities Thus, the PGF of the system state probabilities in the steady-state We need to determine W n , n = 0, 1, ·, c − 1, before P(z) is fully determined. Observe that z = 1 is a trivial root of the denominator of W(z). Using Rouché's theorem we can prove that the denominator has c − 1 other roots inside the open unit ball (i.e., |z| < 1) as shown in the following claim.
. So, we have, for sufficiently small > 0, Consider all the values of z on the contour |z| = 1 + . Using the triangle inequality and (26) we obtain and hence | f (z)| < |g(z) on the contour. Now, since both functions f (z) and g(z) are analytic on the closed disk |z| ≤ 1 + , Rouché's theorem ensures that g(z) and g(z) − f (z) have the same number of zeros in |z| ≤ 1 + , that is, the denominator z c (λ + µ 2 ) − λz c+1 − µ 2 and (λ+µ 2 ) λ z c have the same number of zeros inside the closed disk |z| ≤ 1 + . Letting tend to zero yields the proof of the claim.
These c − 1 roots, ensured by the previous claim, are also the roots of the numerator due to the fact that W(z) is an analytic function on |z| ≤ 1. Replacing these c − 1 roots in the numerator, we obtain c − 1 linear equations with variables W i (i = 0, · · · , c − 1). The c-th linear equation is obtained using the fact P(1) = 1. This equation is equivalent to We now have c linear equations with c unknowns. The linear system of c equations and c variables is solved numerically. Once we have the values of W i , i = 0, · · · , c − 1, the value of W c is obtained from (7).

Measures of Effectiveness
We now calculate some performance measures of the system using the probabilities obtained in this approach. Write W(z) = N(z) D(z) . The expected number of customers in the system in the steady-state is where We may now introduce a cost function to optimize the operations of the system. Let c h be the unit holding cost for each customer in the system, c o be the operating cost per unit of time, c a be the startup cost per unit time for the setup of the server, and c s be the setup cost per busy cycle. Then, the total expected cost per unit of time is where the expected idle period, the expected busy period, and the expected busy cycle are respectively given by Thus,

Illustrative Example
Take for example c = 5, that is, the server adopts the single mode if four customers or less are in the system and the batch mode if five customers or more are in the system. Assume that in the single mode, the service rate is µ 1 = 2, the breakdown rate is α 1 = 0.05, and the repair rate is β 1 = 0.07. In the batch mode, the service rate is µ 2 = 5.5, the breakdown rate is α 2 = 0.08, and the repair rate is β 2 = 0.06. Arrivals at a rate λ = 0.5. The unit costs are c h = 10, c o = 20, c a = 50, and c s = 500. The system of linear equation yields P 0 = 0.7522, P 1 = 0.1875, P 2 = 0.0464, P 3 = 0.0111, P 4 = 0.0023, and P 5 = 0.0002. Also, the average system size is L = 0.3032, and the cost function has value TC = 233.6479.

Numerical Method Using Operators
In this section, we use a different approach to find the probabilities P n . For the sequence of probabilities W n , we define the linear operator D by W n = DW n−1 , ∀n ≥ 1.
Note that composing the operators yields W n+m = D m W n , ∀n, m ≥ 1.
In this case we define, on the interval (0, 1), the function f (z) = µ 2 z c+1 − (λ + µ 2 )z + λ. Similarly to the previous case, the study of the variations and the sign of the values of f gives a unique real root r 3 of f (z) = 0 on (0, 1) which is given explicitly by r 3 := [ λ+µ 2 µ 2 (c+1) ] where d 3 is an arbitrary constant.
Calculation of d 1 , d 2 , and d 3 : Given the values of the roots r 1 , r 2 , and r 3 we are going to determine the values of the constants d 1 , d 2 , and d 3 using the Equation (7), the Equation (9) at n = c, and the summability-to-one condition P(1) = 1. We obtain the following linear system: Calculation of P n : Since P n = W n + F n , we readily have Note that the numerical method gives all the probabilities P n , n ≥ 0, whereas with the analytical method, we only obtain the first c + 1 probabilities P n , 0 ≤ n ≤ c, and the rest of the probabilities P n , n > c, needs to be calculated by successive differentiation.

Measures of Effectiveness
As in the other approach, we calculate the expected number of customers in the system. It is given by The expected idle period is Since we have the explicit form of P 0 , we can find explicitly the mean busy period and the mean busy cycle These measures can be combined to obtain an expression for the total expected cost per unit of time

Illustrative Example
Taking the same parameter as in the previous approach, we find the probabilities P 0 = 0.7496, P 1 = 0.1881, P 2 = 0.0472, P 3 = 0.0118, P 4 = 0.0030, P 5 = 0.0003. We note that the values are remarkably close to the ones obtained in the previous example and the two methods are in total agreement. We also obtain L = 0.3294 and TC = 233.66.

Case Study
The queueing system studied in this paper fits perfectly the following manufacturing situation. Consider a guitar manufacturing factory where guitars can be either handmade or machine made.
The two types of instruments target different market segments. Patrons select their product, basing their choice on different criteria such as quality, value and price, sound, precision, durability, long term repairability, etc. The factory operations manager has implemented a single and batch service strategy as follows: when the number of guitar orders is below some threshold level c, guitars are made by hand, and when it is larger than or equal to c, then guitars are machine made. Machine made guitars are created using a machine to replicate the look and acoustics of an authentic handmade guitar. Because there is minimal labor involved, machine made guitars can be produced quickly, and at a fraction of the price of their handcrafted originals. In this application, guitar orders are the customers. It is assumed that the time between guitar orders is exponential with parameter λ = 4. The server is either the luthier or the machine. Let us assume that the service times and repair times are exponentially distributed.
A handmade guitar will carry a price which reflects its real value in terms of labor and overhead more truly than a factory made one which carries the same price. The former may take two units of time of someone's conscientiously invested time and skill; the latter may take four to seven units of time of intensely repetitive and automated work. Assume the service rate of the luthier is µ −1 1 = 2 units of time and the service rate of the machine is µ −1 2 = 5.5 units of time. Either the luthier or the machine may be unavailable and this happens randomly with respective rates α 1 = 8 and α 2 = 3. Each server becomes available again after a mean time β −1 1 = 1 9 and β −1 2 = 1 6 units of time. The operations manager would like to know the best order level to switch from handmade to machine made guitars. The unit costs of the system are c h = 10, c o = 20, c a = 50, and c s = 500.
Applying the results obtained in the previous section, we calculate the expected total cost TC(c) for successive values of c, starting from c = 1. The variations of the cost are represented in Figure 2. The optimal value of c is found to be c * = 3 and the corresponding optimal cost is TC * = 48.3038. In managerial terms for the operations manager, this means that the optimal policy is to have the luthier make the guitars by hand as long as there are less than three orders in line. If this number is grater than or equal to 3, then guitars should be made using the machine. In case there is some uncertainty about some of the parameters, a sensitivity analysis can be conducted to assess the effect of this uncertainty on the optimal policy. For example, suppose there is an uncertainty about α 1 , the failure rate when the server is in the working state. Then optimal measures can be calculated for various values of α 1 . A sample calculations is shown in Table 1 where we calculate the probability of no orders P 0 , the average number of orders L, and the optimal expected cost per unit of time TC * . The effect of other parameters can be assessed in a similar way.

Conclusions
The Markovian queueing system considered in this paper is characterized by a flexible server that adapts to the queue length by switching from a single service to a bulk service when the queue length is too large and from bulk service to single service when the queue length is too small. The server is unreliable and may break down while providing service. Different parameters depend on the service discipline applied. We calculated the system size steady-state probabilities in terms of their probability generating function and using linear operators. The two methods comply with each other. An application to a case study is also provided.
There are various ways this work can be further developed. For example, bulk arrivals instead of single arrivals could be examined. Also general distributions could be assumed for the various processes considered. Server vacations, either working or not, and various threshold policies such as N-, T-, or D-policies could also be taken into account.