Abstract
Consider a closed cyclic queueing model that consists of two nodes and a total of customers. Each node buffer can accommodate all customers. Node 1 has servers, each having an exponential service time with rate . The second node consists of a single server with a general service time distribution function . The well-known machine repair model with spares, where a set of identical machines, , is served by a single repair person, is a key application of this model. This model has several other applications in performance evaluation, manufacturing, computer networks, and in reliability studies as it can be easily used to compute system availability. In this article, we give an efficient algorithm to derive an exact solution for the steady state system size probabilities. Our approach is based on developing an efficient polynomial convolution method to compute the transition probabilities of a birth process over node 2 service times and solving an imbedded Markov chain at node 2 service completion epochs. This is a significant improvement over the exponential algorithm given in an earlier paper. Numerical examples are given to demonstrate the performance of our method.
1. Introduction
Computing the steady-state probability distribution of a non-Markovian two-node closed queueing cyclic network is known to be computationally challenging. In this article, we propose a new convolution method that is significantly more efficient than existing algorithms. Specifically, we consider a two-node closed cyclic queueing network with a total of customers. The first node consists of parallel identical servers having independent and identically distributed (i.i.d.) exponential service times with rate . The second node consists of a single server with i.i.d. service requirements , , having a general distribution function . Waiting space at each node is sufficient to accommodate all customers and the service discipline is FCFS. Figure 1 shows the diagram of the model.
A key contribution of this paper is to develop an easy to understand efficient polynomial algorithm to compute the stationary distribution of the number of customers at node 2 and to determine the measures of performance of this two-node network model. For this purpose, we develop and solve an imbedded Markov Chain (MC) with states given by the number of customers in the second node immediately after departure epochs. Such an approach is standard in analyzing the general single server model with Markovian arrivals, . In our model, this is equivalent to the special case of . The fact that can be larger than complicates the development of the imbedded MC in our model since the transition rate to the second node is not as smooth as in the case of and . Except for Maddah and El-Taha [1], no other article utilizes an imbedded Markov chain approach for the closed cyclic two-node network model as the one we consider.
To the best of our knowledge, only a few articles have addressed this problem before. Gupta and Rao [2] applied the supplementary variable technique to develop the system size probabilities for a model similar to ours. Gupta and Rao [2] solution approach draws heavily on transform techniques and difference equations and requires higher order derivatives of the Laplace transform of the node 2 service time distribution. Recently, Oz et al. [3] utilized conditional residual service times and a recursive scheme to compute the stationary distribution function. Maddah and El-Taha [1] use the imbedded Markov chain approach representing node 2 system size probabilities at departure epochs. They compute the transition probabilities using an exponential algorithm. This paper’s contribution is a significant improvement over the results of Maddah and El-Taha [1] by giving a more efficient polynomial algorithm to compute transition probabilities. The imbedded Markov chain technique has been successfully applied to solve the related machine repairman model with no spares, (see, for example, Van Hoorn [4], Saaty [5], and Takagi [6]). The papers by Boxma [7] and Daduna [8] deal with two-node network with single server at each node. Boxma and Daduna [9] discuss sojourn time distributions and approximate methods for two-node queueing networks. A more detailed literature review can be found in Maddah and El-Taha [1].
Our model applies to the machine repairmen model with spare machines under the following assumptions. Machine times to failure are iid exponential random variables with mean . The repair facility is staffed by a single server with iid repair requirements , having a general distribution function . The waiting room at both the repair facility and the production floor can accommodate all machines (no blocking), and the service discipline is FCFS.
The rest of the article is organized as follows. In Section 2, we introduce a birth process and evaluate its transition probabilities. In Section 3, we determine the one-step transition probabilities of the Markov chain imbedded at the departure instants. In Section 4, we solve for the system size probabilities at node 2 departure instants, the time average stationary distribution, and measures of performance of the model. Specifically, we give a stable algorithm that solves for the stationary distribution of the imbedded Markov chain and the time-average probabilities. In Section 5, we give small size numerical examples and solve large-size problems. In Appendix A we give a second direct proof of our main result in Theorem 3, and in Appendix B, we describe in detail the algorithm to compute the stationary probabilities.
2. Birth Process
Consider the two-node model with . Let be a pure birth process representing the number of births (node 1 service completions) during one node 2 service time of length . Birth rates (i.e., arrival rates at node 2) depend on the number of customers at node 2 and are given by , where is the number of customers at node 2 immediately after a departure. That is a pure birth process is due to the fact that the time between two node 2 consecutive arrivals (machine failures) is exponential, possibly with different rates. Let be the transition probabilities of the birth process defined as where represents the number of node 2 customers at time zero.
Let be a sequence of independent random variables with distribution functions . Here, represents interevent times. Let be the d.f. of which is the -fold convolution of .
Lemma 1. Let be a simple point process defined as , such that and as . Moreover, let the times between events be independent but not necessarily identically distributed. Then, .
Proof (Similar to the renewal case). The difference between Lemma 1 and the well-known renewal case is that we relax the assumption that interevent times are identically distributed. Noting that , we obtain This completes the proof of the lemma.☐
We will also need the following convolution result. Let , where and are nonnegative random variables. Then,
Now, we state our main result.
Theorem 2. For ; , where , over an empty set is zero, and over an empty set is one.
Proof. Let , where ; that is, is an exponential random variable with parameter , so that . Therefore, the pdf of is given by Moreover, let , where Thus, hypoexponential (see Ross [10], p 299), and where , and where a product over an empty set is one. Therefore, Now, by Lemma 1 and Equation (3) Therefore, Noting that for , we obtain Therefore, Note that the second term inside the brackets of the equation is the same as the first when . This completes the proof.☐
For the special case , , we obtain which is consistent with the approach given by Maddah and El-Taha [1].
By referring to Figure 2, we note that in regions I and III are, respectively, similar to the and the finite population models. See Maddah and El-Taha [1] Lemma 2.2(i) and Lemma 2.2 for derivations of the in these two regions. The in Region II, given by Theorem 2, is the most challenging since the transition from to is not smooth as in the other regions. Region IV is where is .
2.1. Remarks on Complexity
We note here that the advantage of Theorem 4.4 of Maddah and El-Taha [1] is that for the first time it gives a closed form expression, based on a multidimensional triangular geometric sum (MTGS) result, for computing . That was an improvement over the supplementary approach that did not have a closed form expression and required the computation of higher order derivatives of the Laplace transform of the nonexponential service times at node 2. Our convolution method here is a significant improvement on Maddah and El-Taha [1].
We compare the complexity of computing using the convolution method given in Theorem 2 and the method based on the MTGS in Theorem 4.4 of Maddah and El-Taha [1]. Now, the complexity of the computations of the convolution method in 2.2 is on the order of , while the complexity of the computations of the MTGS method in Theorem 4.4 in Maddah and El-Taha [1] is on the order of. For example, for a problem with and , the complexity of the convolution method is on the order of computations while that of MTGS method is of the order of computations. Moreover, for a problem with and , the convolution method has a complexity of a manageable order of computations, while that of the MTGS has the unmanageable order of computations.
3. One-Step Transition Probabilities
In this section, we determine the transition probabilities, the system size probabilities, and measures of performance of the two-node cyclic model described above. We develop an embedded Markov chain and then solve for the system size departure epochs’ probabilities. Let be a stochastic process that represents the number of customers present at node 2 immediately following a departure (i.e., be the number of machines in the repair facility upon the repair completion of the machine). Then, is related to as follows: where is the number of failures during the repair time of the machine. Define the departure instants distribution as
The transition probabilities of , defined as , are given by
and for , for all ; otherwise, .
Let ; then , and where are the transition probabilities of . Now, the transition probability matrix of the Markov chain is evaluated utilizing Equation (18) above.
Let be the Laplace-Stieltjes transform (LST) of the service-times distribution function . Moreover, let , where is the derivative of , with . Then, . Now, we have the following result.
Theorem 3. The one-step transition probabilities are given by (i)for , (ii)for , (iii)for , where, .
Proof. The proofs of (i) and (ii) follow from Lemma 2.2(i) and Lemma 2.2 of Maddah and El-Taha [1]. Because , we have Use (18) and simplify to obtain Now, for all , write where the over an empty an set is zero, and over an empty set is one. This completes the proof.☐
In contrast to , the term is always positive.
Remark 4. For computational considerations and to avoid the propagation of error due to subtractions, we write (20) as where is the notation for the greatest integer. Similarly, we write (21) as
Although these expressions appear more complicated, they are more computationally efficient than the original ones. In expressions (25) and (26), we subtract only once at the end versus alternating between addition and subtraction in (20) and (21).
Computing the one-step transition probabilities requires the evaluation of the derivatives of the LST of the service time distribution functions. This is explained further in Section 5 where we give closed form expressions for the derivatives. Equation (20) allows the evaluation of with ease for service time distributions with closed form Laplace transforms. The second expression of in Theorem 3 is similar to the case of the model for which some efficient computation techniques have been established. The third expression in (21) is obtained using Theorem 2.
4. Stationary Distributions
In this section, we solve the imbedded Markov chain using a stable efficient approach to evaluate the departure-time stationary distribution function, then establish a relationship between the departure-time probabilities and time-average probabilities which we use to compute the time-average distribution.
4.1. Departure-Time Probabilities
Having determined the transition probability matrix for the Markov chain , the limiting probabilities can be evaluated by solving the global balance equations . Instead, we give a more stable recursive algorithm to solve for the departure-time probabilities . Using the fact that probability flow across cuts balances (Kelly [11], Lemma 1.4), one can show that for , which implies where .
Now, can be computed recursively using (28), which is more numerically stable than the method based on the global balance equations as shown by Stidham [12].
4.2. Time-Average Probabilities
In this subsection, we are interested in node 2 time-average stationary distribution defined as
Moreover, we are interested in determining the performance measures of the system.
Let be node 2 mean service time, for , and otherwise. Then, the throughput is given by
The system size probabilities, , are given by and .
Moreover, the mean number of node 2 customers is given by Little’s law implies that node 2 waiting time is . Furthermore, the mean delay is and the mean number of node 1 customers is . In the next section, we compute the stationary distribution and calculate , , and using several service time distributions.
5. Applications
In this section, we consider both small-scale and large-scale examples. For both small-scale and large-scale applications, we focus on service time distributions whose LST have closed form multiple derivatives. We choose four distribution functions with coefficients of variation that vary from to infinity. Specifically, we select the deterministic, Erlang, exponential, and the hyperexponential distribution functions. Let be the of the service times. Recall that we require the service time distribution function to have a mean . Here, we give explicit forms for the derivatives of the LST of the service time distribution functions.
Deterministic. In this case, we assume that w.p. 1, so that
Exponential. Here , so that
Erlang. The density function for . Here,
Hyperexponential. The two-phase hyperexponential is a mixture of two exponential distribution functions and has a density function of the form . Therefore,
For the above distributions, we recursively compute the stationary distributions and for the examples in the next subsection.
5.1. Small-Scale Applications
Here, we present a small-scale, two-node network with a total of customers and servers at the first node, so that . Now, referring to Theorem 3, we write the one-step transition probabilities in the more simplified convenient form as follows. (i)For Theorem 3 (i), where , (i.e., and ), we use (ii)For Theorem 3 (ii), where , (i.e., and ), we use (iii)For Theorem 3 (iii), where (i.e., and ), we use the convolution result where
Furthermore, for and , we simply use . Moreover, where .
When , we use for all . Intuitively, if is the number of customers at node 2 immediately after a departure, then the transition probabilities when must be the same as the probabilities when , because in the later case, node 2 becomes empty, and in the former case, node 2 is already empty. Therefore, the one-step transition matrix is given by
Next, we specialize the same example above using deterministic service times.
5.1.1. Deterministic Service Times
Let node 2 service times be a deterministic probability density function, , and otherwise (i.e., ); thus, we have therefore, and This gives the following: and Substitute in the general service time matrix to get the corresponding one-step transition matrix for deterministic service times.
As an illustration, consider the case of a machine repair model with deterministic repair times; that is, the repair time is equal to . The system owns six machines, with a demand of three machines, so that and , so that . Thus, our transition matrix is as follows.
Similar to the deterministic distribution function, we can generate one-step transition matrices for other distributions including the exponential, Erlang, and hyperexponential.
5.1.2. A Reliability Example: System Availability
The system consists of components subject to failure and repair. In this example, we consider the availability of a -out-of- system with spares, exponential failure, and general repair times. In this system, at least out of the components must function for the system to function. The key measure of performance here is to obtain the system availability which is . The one-step transition probabilities for this system are computed in the same way as our model except that here the state space is restricted to . This reliability model is considered in example 4 by Maddah and El-Tha [1]. With our approach, we can solve larger reliability problems.
5.2. Numerical Results
Here, we provide numerical results for small- and large-scale examples. We start with small-scale applications. Consider the following service time node 2 distributions: deterministic, exponential, Erlang, and two-phase hyperexponential () distribution discussed earlier in this section. We fix the mean of the distribution to , i.e., , for all four distributions. For , we choose to guarantee high variability. Let be the offered load at the repair server. Applying our algorithm to the four distributions for the system with customers, , and , we obtain the results in Tables 1 and 2.
To verify our numerical computations, results with exponential service/repair times were compared with similar values obtained from the classical approach based on a birth death process. Results from both approaches were the same for all cases as expected. For other distributions, we utilized simulation to verify the results. All results were valid when compared with simulation outputs with an average run length of one million repair completions.
For large-size problems, the algorithm will become unstable at some point. We tested the algorithm for a large number of cases with various parameters. We report on two cases where the system is stable. The first example is for , , so that . We set the service rate at node 2 to . Service rate at node 1 (or machine failure rate) is set to so that we have a reasonable system throughput or utilization of node 2 server. Additionally, note that the offered load per active node 1 server is . All steps of the algorithm are the same for all repair-time distribution functions. We verified our results by simulation where we simulated 10 million node 2 service completions. As expected, the results in Table 3 show that the system has the highest throughput for the deterministic and the lowest throughput for the hyperexponential service time distribution functions. Note that , is node 2 effective traffic intensity.
For the second example, we use , , so that . We set the service rate at node 2 to . Service rate at node 1 (or machine failure rate) is set to so that we have a reasonable system throughput or utilization of node 2 server. Furthermore, note that the offered load per active node 1 server is . A summary of the results is given in Table 4.
Note that changing from in Table 3 to in Table 4 (thus decreasing variability) for the hyperexponential service times reduces significantly. In Figures 3 and 4, we present graphically the and the of the number of customers at node 2 using the four service time distribution functions. As expected, the hyperexponential has the heaviest tail.
The hyperexponential presents the biggest challenge. To verify the correctness of our computations, we compared the of the number in the system using both the algorithm and the simulation output. We had a perfect match.
In conclusion, our algorithm represents a significant improvement over Maddah and El-Tha [1] by solving accurately much larger problems.
Appendix
A. An additional Proof
A.1. A Second Direct Proof of Theorem 3 (iii)
Recall that is the number of customers at node 2 immediately after a service completion. Given that , define as the time required for arrivals, which takes the state from to , and arrivals for the remaining time.
Now, let be the time it takes to go from to , and . Now, is an exponential pdf with parameter , and is the sum of random variables. Therefore, has an Erlang distribution function .
That is, given , the time for arrivals follows an Erlang distribution function with pdf
Now, going from to is given by
This leads to the equation
Note that which leads to
where we replace by in Step 4. Now, noting that is a gamma density function, we obtain (e.g., Tijms [13], page 442)
Thus,
Similar to the proof of Theorem 3 (iii)
Refer to Theorem 3 (iii) and expand as a factorial.
Substitute in (21) to obtain
which is the same as (A.8). This provides an alternative path to prove Theorem 3 (iii).
B. Algorithm
The convolution algorithm presented here is used to compute the stationary probabilities and performance measures for both the small-scale and large-scale examples. Our computational methodology follows that is provided in Sections 3 and 4 and the first part of Section 5. We point out that the solution method is simple and consists of a few steps: (i)Compute the one-step transition probabilities ’s.(ii)Compute recursively using Equation (28).(iii)Compute using Equation (31).(iv)Compute the desired performance measures.
Here, we give a detailed algorithm that provides a blueprint for a reader interested in coding this method using any programming language.
Algorithm B.1. Initialization. Let be the total number of customers in the system, be the number of parallel servers at node 1, be the service rate at node 1, and be the mean service time at node 2. Moreover, for each specified service time distribution function, let be such that
Deterministic:
Exponential:
Erlang (two-phase):
Hyperexponential:
Note that for , .
(1)For each specified service time distribution, compute such that
(i)for , compute ;(ii)for , compute ; and(iii)for , compute , and compute .Computations in (i), (ii), and (iii) are to be used with (4), (5), and (6). below respectively.
(2)Compute for , using
where over an empty set is zero, and over an empty set is one.(3)Set for all , .(4)Compute for , using
(5)Compute for , using
where is the greatest integer in .(6)Compute for , using
(7)Compute for ,(8)Compute , for ,(9)Define for , (10)Compute for , using
(11)Set (12)Compute for recursively using
(13)Compute for by normalizing using
where
(14)Compute , and .(15)Compute for using
and
Compute performance measures using and .
Data Availability
The data used in the article is generated by a python program that can be made available upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.