An Efficient Convolution Algorithm for the Non-Markovian Two-Node Cyclic Network

Consider a closed cyclic queueing model that consists of two nodes and a total ofM customers. Each node buffer can accommodate all M customers. Node 1 has N ≤M servers, each having an exponential service time with rate λ. The second node consists of a single server with a general service time distribution function Bð:Þ. The well-known machine repair model with spares, where a set of identical machines, N , 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.


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 M customers. The first node consists of N ≤ M 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 S i , i = 1, 2, ⋯, having a general distribution function Bð:Þ. Waiting space at each node is sufficient to accommodate all M 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, M/G/1. In our model, this is equivalent to the special case of M = N. The fact that M can be larger than N 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 M/G/1 and M/G/1//N. 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, M/G/1//N (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 1/λ. The repair facility is staffed by a single server with iid repair requirements S i , i = 1, 2, ⋯ having a general distribution function Bð:Þ. The waiting room at both the repair facility and the production floor can accommodate all M 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.

Birth Process
Consider the two-node model with M = N + Y. Let fYðtÞ, t ≥ 0g be a pure birth process representing the number of births (node 1 service completions) during one node 2 service time of length t. Birth rates (i.e., arrival rates at node 2) depend on the number of customers at node 2 and are given by λ i = min fN, ðM − iÞgλ, where i is the number of customers at node 2 immediately after a departure. That fYðtÞ, t ≥ 0g 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 Q ij ðtÞ ; j ≥ i be the transition probabilities of the birth process fYðtÞ, t ≥ 0g defined as where Yð0Þ = i represents the number of node 2 customers at time zero. Let fX i , i ≥ 1g be a sequence of independent random variables with distribution functions F i . Here, X i represents interevent times. Let F n ðtÞ be the d.f. of S n = X 1 + ⋯ + X n which is the n-fold convolution of F 1 , ⋯, F n . Lemma 1. Let fNðtÞ, t ≥ 0g be a simple point process defined as NðtÞ = max fn : S n ≤ tg, such that Nð0Þ = 0 and NðtÞ ⟶ ∞ as t ⟶ ∞. Moreover, let the times between events be independent but not necessarily identically distributed. Then, PðNðtÞ = nÞ = F n ðtÞ − F n+1 ðtÞ = F n+1 ðtÞ − F n ðtÞ.
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 fNðtÞ ≥ ng ⇔ fS n ≤ tg, we obtain This completes the proof of the lemma.☐ We will also need the following convolution result. Let Y = X + Z, where X and Z are nonnegative random variables. Then, Now, we state our main result.
where C k,j−Y = Q j−Y m=1 m≠k ððN − mÞ/ðk − mÞÞ, ∑ over an empty set is zero, and Q over an empty set is one.
Proof. LetX i = X i + ⋯ + X Y , where X k~e xp ðNλÞ, k = i, ⋯, Y; that is, X k is an exponential random variable with  Journal of Applied Mathematics parameter Nλ, so thatX i~E rlangðY − i + 1, NλÞ. Therefore, the pdf ofX i is given by Moreover, letZ j = Z 1 + ⋯Z j , where Z k~e xp ððN − kÞλÞ, k = 1, ⋯, j − Y: Thus,Z j~h ypoexponential (see Ross [10], p 299), and where λ k = ðN − kÞλ, and where a product over an empty set is one. Therefore, Now, by Lemma 1 and Equation (3) Therefore, Noting that for k < j − Y, we obtain Therefore, Note that the second term inside the brackets of the equation is the same as the first when k = j − Y. This completes the proof.☐ For the special case j = Y + 1, Y ≥ 2, we obtain 3 Journal of Applied Mathematics which is consistent with the approach given by Maddah and El-Taha [1].
By referring to Figure 2, we note that Q ij ðtÞ in regions I and III are, respectively, similar to the M/G/1 and the M/G/ 1/N finite population models. See Maddah and El-Taha [1] Lemma 2.2(i) and Lemma 2.2 for derivations of the Q ij ðtÞ in these two regions. The Q ij ðtÞ in Region II, given by Theorem 2, is the most challenging since the transition from i to j is not smooth as in the other regions. Region IV is where Q ij ðtÞ is 0.

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 Q ij ðtÞ. 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 Q ij ðtÞ using the convolution method given in Theorem 2 and the method based on the MTGS in Theorem 4.4 of Maddah and El-Taha

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 fX n , n ≥ 1g 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 n th machine). Then, X n is related to X n−1 as follows: where A n is the number of failures during the repair time of the n th machine. Define the departure instants distribution The transition probabilities of fX n , n ≥ 1g, defined as pði, and for i = 0, pð0, jÞ = pð1, jÞ for all j; otherwise, pði, jÞ = 0, j < i − 1.

Journal of Applied Mathematics
Proof. The proofs of (i) and (ii) follow from Lemma 2.2(i) and Lemma 2.2 of Maddah and El-Taha [1]. Because pði, j ; tÞ = Q i,j+1 ðtÞ, we have Use (18) and simplify to obtain Now, for all where the ∑ over an empty an set is zero, and Q over an empty set is one. This completes the proof.☐ In contrast to C k,j−Y+1 , the term C a k,j−Y+1 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

Journal of Applied Mathematics
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 pði, jÞ 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

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.

Departure-Time Probabilities.
Having determined the transition probability matrix for the Markov chain fX n , n = 0, ⋯g, the limiting probabilities πðiÞ, i = 0, ⋯, M − 1 can be evaluated by solving the global balance equations πP = π ; ∑ M−1 i=0 πðiÞ = 1. 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 j = 1, 2, ⋯, M − 1, which implies where aðk, jÞ = ∑ M−1 i=j pðk, iÞ. Now, fπð:Þg 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].

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 E½S = 1/μ be node 2 mean service time, λ i = Nλ for i ≤ Y, and λ i = ðM − iÞλ otherwise. Then, the throughput λ f is given by   Journal of Applied Mathematics The system size probabilities, pðiÞ, i = 0, ⋯, M, are given by and pðMÞ = 1 − ∑ M−1 i=0 pðiÞ. Moreover, the mean number of node 2 customers L is given by L = ∑ M i=1 ipðiÞ: Little's law implies that node 2 waiting time is W = L/λ f . Furthermore, the mean delay is W q = W − E½S, and the mean number of node 1 customers is M − L. In the next section, we compute the stationary distribution and calculate λ f , L, and W using several service time distributions.

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 0 to infinity. Specifically, we select the deterministic, Erlang, exponential, and the hyperexponential distribution functions. Let bðtÞ be the p:d:f : of the service times. Recall that we require the service time distribution function to have a mean E½B = 1/μ. Here, we give explicit forms for the derivatives of the LST of the service time distribution functions.

Small-Scale Applications.
Here, we present a small-scale, two-node network with a total of M = 6 customers and N = 3 servers at the first node, so that Y = 3. Now, referring to Theorem 3, we write the one-step transition probabilities in the more simplified convenient form as follows.
As an illustration, consider the case of a machine repair model with deterministic repair times; that is, the repair time is equal to b w:p:1. The system owns six machines, with a demand of three machines, so that M = 6 and N = 3, so that Y = 3. Thus, our transition matrix is as follows.

A Reliability
Example: System Availability. The system consists of M components subject to failure and repair. In this example, we consider the availability of a K-out-of-M system with Y spares, exponential failure, and general repair times. In this system, at least K out of the M components must function for the system to function. The key measure of performance here is to obtain the system availability which is 1 − PðM − K + 1Þ. 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 f0, 1,⋯,M − K + 1g. This reliability model is considered in example 4 by Maddah and El-Tha [1]. With our approach, we can solve larger reliability problems.

Numerical Results.
Here, we provide numerical results for small-and large-scale examples. We start with smallscale applications. Consider the following service time node 2 distributions: deterministic, exponential, Erlang, and twophase hyperexponential (H 2 ) distribution discussed earlier in this section. We fix the mean of the distribution to 1:0, i.e., μ = 1, for all four distributions. For H 2 , we choose p = 0:005 to guarantee high variability. Let ρ = λ/μ be the offered load at the repair server. Applying our algorithm to the four distributions for the system with M = 6 customers, N = 3, and Y = 3, 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 M = 100, N = 80, so that Y = M − N = 20. We set the service rate at node 2 to μ = 1. Service rate at node 1 (or machine failure rate) is set to λ = :012 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 ρ = λ/μ = :012. 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 λ f for the deterministic and the lowest throughput for the hyperexponential service time distribution functions. Note that λ f = ρ e , is node 2 effective traffic intensity.
For the second example, we use M = 60, N = 40, so that Y = 20. We set the service rate at node 2 to μ = 1. Service rate at node 1 (or machine failure rate) is set to λ = :025 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 ρ = λ/μ = :025. A summary of the results is given in Table 4.
Note that changing p from :005 in Table 3 to :125 in Table 4 (thus decreasing variability) for the hyperexponential service times reduces L significantly. In Figures 3 and 4, we present graphically the pmf and the cdf 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 cdf 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.