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

Authors: Muhammad El-Taha; Thomas Michaud
Journal: Journal of Applied Mathematics
 (2021)

---

## Abstract

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.


---

## Body

## 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 ofM 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 Si, 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.Figure 1
Finite population 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 areiid exponential random variables with mean 1/λ. The repair facility is staffed by a single server with iid repair requirements Si, 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 Section2, 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 withM=N+Y. Let Yt,t≥0 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=minN,M−iλ, where i is the number of customers at node 2 immediately after a departure. That Yt,t≥0 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 Qijt;j≥i be the transition probabilities of the birth process Yt,t≥0 defined as
(1)Qijt=PYt=j−iY0=i,where Y0=i represents the number of node 2 customers at time zero.LetXi,i≥1 be a sequence of independent random variables with distribution functions Fi. Here, Xi represents interevent times. Let Fnt be the d.f. of Sn=X1+⋯+Xn which is the n-fold convolution of F1,⋯,Fn.Lemma 1.
LetNt,t≥0 be a simple point process defined as Nt=maxn:Sn≤t, such that N0=0 and Nt⟶∞ as t⟶∞. Moreover, let the times between events be independent but not necessarily identically distributed. Then, PNt=n=Fnt−Fn+1t=F¯n+1t−F¯nt.Proof (Similar to the renewal case).
The difference between Lemma1 and the well-known renewal case is that we relax the assumption that interevent times are identically distributed. Noting that Nt≥n⇔Sn≤t, we obtain
(2)PNt=n=PNt≥n−PNt≥n+1=PSn≤t−PSn+1≤t=Fnt−Fn+1t.
This completes the proof of the lemma.☐We will also need the following convolution result. LetY=X+Z, where X and Z are nonnegative random variables. Then,
(3)F¯Yt=∫0tfXsF¯Zt−sds.Now, we state our main result.Theorem 2.
Fori=1,⋯Y−1; j=Y+1,⋯M−1,
(4)Qijt=NY−i+1e−Nλt∑k=1j−YN−kCk,j−YN−j+YkY−i+1×ekλt−∑r=0Y−ikλtrr!,where Ck,j−Y=∏m=1m≠kj−YN−m/k−m, ∑ over an empty set is zero, and ∏ over an empty set is one.Proof.
LetX~i=Xi+⋯+XY, where Xk~expNλ,k=i,⋯,Y; that is, Xk is an exponential random variable with parameter Nλ, so that X~i~ErlangY−i+1,Nλ. Therefore, the pdf of X~i is given by
(5)fX~is=NλsY−iNλe−NλsY−i!.
Moreover, letZ~j=Z1+⋯Zj, where Zk~expN−kλ,k=1,⋯,j−Y. Thus, Z~j~ hypoexponential (see Ross [10], p 299), and
(6)F¯Z~jt=∑k=1j−YCk,j−Ye−λkt,where λk=N−kλ, and
(7)Ck,j−Y=∏m=1m≠kj−Yλmλm−λk=∏m=1m≠kj−YN−mλN−mλ−N−kλ=∏m=1m≠kj−YN−mk−m,where a product over an empty set is one. Therefore,
(8)F¯X~i+Z~jt=∫0tNλsY−iNλe−NλsY−i!×∑k=1j−YCk,j−Ye−N−kλt−sds=∑k=1j−YCk,j−YNλY−i+1×e−N−kλt∫0tsY−ie−NλseN−kλsY−i!ds=∑k=1j−YCk,j−YNλY−i+1×e−N−kλt∫0tsY−iY−i!e−kλsds=∑k=1j−YCk,j−YNλY−i+1e−N−kλtkλY−i+1∫0tkλkλsY−iY−i!e−kλsds=∑k=1j−YCk,j−YNλY−i+1e−N−kλtkλY−i+11−∑r=0Y−ikλtre−kλtr!=NY−i+1e−Nλt∑k=1j−YCk,j−YekλtkY−i+1×1−e−kλt∑r=0Y−ikλtrr!=NY−i+1e−Nλt∑k=1j−YCk,j−YkY−i+1ekλt−∑r=0Y−ikλtrr!.
Now, by Lemma1 and Equation (3)
(9)Qijt=PNt=j−i∣N0=i=F¯X~i+Z~jt−F¯X~i+Z~j−1t.
Therefore,(10)Qijt=NY−i+1e−Nλt∑k=1j−YCk,j−YkY−i+1×ekλt−∑r=0Y−ikλtrr!−∑k=1j−Y−1Ck,j−Y−1kY−i+1ekλt−∑r=0Y−ikλtrr!=NY−i+1e−Nλt∑k=1j−Y−1Ck,j−YkY−i+1ekλt−∑r=0Y−ikλtrr!+Cj−Y,j−Yj−YY−i+1ej−Yλt−∑r=0Y−ij−Yλtrr!−∑k=1j−Y−1Ck,j−Y−1kY−i+1ekλt−∑r=0Y−ikλtrr!=NY−i+1e−Nλt∑k=1j−Y−1Ck,j−Y−Ck,j−Y−1kY−i+1ekλt−∑r=0Y−ikλtrr!+Cj−Y,j−Yj−YY−i−1ej−Yλt−∑r=0Y−ij−Yλtrr!.
Noting that fork<j−Y,
(11)Ck,j−Y=Ck,j−Y−1λj−Yλj−Y−λk,we obtain
(12)Ck,j−Y−Ck,j−Y−1=Ck,j−Y1−λj−Y−λkλj−Y=Ck,j−Yλkλj−Y=Ck,j−YN−kN−j+Y.
Therefore,(13)Qijt=NY−i+1e−Nλt∑k=1j−Y−1N−kCk,j−YN−j+YkY−i+1×ekλt−∑r=0Y−ikλtrr!+Cj−Y,j−Yj−YY−i+1ej−Yλt−∑r=0Y−ij−Yλtrr!.
Note that the second term inside the brackets of the equation is the same as the first whenk=j−Y. This completes the proof.☐For the special casej=Y+1, Y≥2, we obtain
(14)Qi,Y+1t=NY−i+1eNλteλt−∑r=0Y−iλtrr!,which is consistent with the approach given by Maddah and El-Taha [1].By referring to Figure2, we note that Qijt 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 Qijt in these two regions. The Qijt 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 Qijt is 0.Figure 2
Birth process transition probabilities.


### 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 Qijt. 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 computingQijt 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 ONYY+1/2, while the complexity of the computations of the MTGS method in Theorem 4.4 in Maddah and El-Taha [1] is on the order ofOM−2Y−1. For example, for a problem with N=20 and Y=10, the complexity of the convolution method is on the order of 1,100 computations while that of MTGS method is of the order of 6,906,900 computations. Moreover, for a problem with N=50 and Y=50, the convolution method has a complexity of a manageable order of 61,250 computations, while that of the MTGS has the unmanageable order of 9.4×10153 computations.


## 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 Qijt. 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 computingQijt 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 ONYY+1/2, while the complexity of the computations of the MTGS method in Theorem 4.4 in Maddah and El-Taha [1] is on the order ofOM−2Y−1. For example, for a problem with N=20 and Y=10, the complexity of the convolution method is on the order of 1,100 computations while that of MTGS method is of the order of 6,906,900 computations. Moreover, for a problem with N=50 and Y=50, the convolution method has a complexity of a manageable order of 61,250 computations, while that of the MTGS has the unmanageable order of 9.4×10153 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. LetXn,n≥1 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 nth machine). Then, Xn is related to Xn−1 as follows:
(15)Xn=Xn−1+An−1,ifXn−1>0,An,ifXn−1=0,where An is the number of failures during the repair time of the nth machine. Define the departure instants distribution πi,i=0,⋯,M−1 as
(16)πi=limn⟶∞PXn=i.The transition probabilities ofXn,n≥1, defined as pi,j=PXn=j∣Xn−1=i,i=0,⋯,M−1;j=0,⋯,M−1, are given by
(17)pi,j=PAn=j−i+1∣Xn−1=i,and fori=0, p0,j=p1,j for all j; otherwise, pi,j=0,j<i−1.Letpi,j;t=PYt=j+1∣Y0=i; then pi,j;t=Qi,j+1t, and
(18)pi,j=∫0∞PYt=j+1∣Y0=idBt=∫0∞Qi,j+1tdBt,where Qijt are the transition probabilities of Yt,t≥0. Now, the transition probability matrix P=pi,j of the Markov chain Xn,n=0,1,⋯,M−1 is evaluated utilizing Equation (18) above.LetB∗s=∫0∞e−stdBt be the Laplace-Stieltjes transform (LST) of the service-times distribution function Bt. Moreover, let Bn∗s=−1ndnB∗s/dsn, where dnB∗s/dsn is the nth derivative of B∗s, n=0,1,⋯, with B0∗s=B∗s. Then, Bn∗s=∫0∞tne−stdBt. Now, we have the following result.Theorem 3.
The one-step transition probabilities are given by(i)
fori=1,⋯,Y−1;j=i−1,⋯Y−1,
(19)pi,j=Nλj−i+1Bj−i+1∗Nλj−i+1!,(ii)
fori=Y,⋯,M−1;j=i−1,⋯,M−2,
(20)pi,j=M−ij−i+1∑r=0j−i+1−1rj−i+1rB∗M−j−1+rλ,(iii)
fori=1,⋯,Y−1;j=Y,⋯,M−2,
(21)pi,j=NY−i+1∑k=1j−Y+1−1j−Y−k+1N−kCk,j−Y+1aN+Y−j−1kY−i+1B∗N−kλ−∑r=0Y−ikλrBr∗Nλr!,where, Ck,j−Y+1a=∏m=1k−1N−m/k−m×∏m=k+1j−Y+1N−m/k−m.Proof.
The proofs of (i) and (ii) follow from Lemma 2.2(i) and Lemma 2.2 of Maddah and El-Taha [1]. Because pi,j;t=Qi,j+1t, we have
(22)pi,j;t=NY−i+1e−Nλt∑k=1j−Y+1N−kCk,j−Y+1N−j+Y+1kY−i+1×ekλt−∑r=0Y−ikλtrr!.
Use (18) and simplify to obtain
(23)pi,j=NY−i+1∑k=1j−Y+1N−kCk,j−Y+1N+Y−j−1kY−i+1B∗N−kλ−∑r=0Y−ikλrBr∗Nλr!.
Now, for all1≤k≤j−Y+1, write
(24)Ck,j−Y+1=∏m=1k−1N−mk−m×∏m=k+1j−Y+1N−mk−m=−1j−Y−k+1∏m=1k−1N−mk−m×∏m=k+1j−Y+1N−mm−k=−1j−Y−k+1Ck,j−Y+1a,where the ∑ over an empty an set is zero, and ∏ over an empty set is one. This completes the proof.☐In contrast toCk,j−Y+1, the term Ck,j−Y+1a is always positive.Remark 4.
For computational considerations and to avoid the propagation of error due to subtractions, we write (20) as
(25)pi,j=M−i!M−j−1!×∑u=0j−i+1/2B∗M−j−1+2uλ2u!j−i+1−2u!−∑u=0j−i/2B∗M−j+2uλ2u+1!j−i−2u!,where . is the notation for the greatest integer. Similarly, we write (21) as
(26)pi,j=NY−i+1N+Y−j−1×∑u=1j−Y+2/2N−2u+1C2u−1,j−Y+1a2u−1Y−i+1B∗N−2u+1λ−∑r=0Y−i2u−1λrBr∗Nλr!−∑u=1j−Y+1/2N−2uC2u,j−Y+1a2uY−i+1B∗N−2uλ−∑r=0Y−i2uλrBr∗Nλr!.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 probabilitiespi,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 (20) allows the evaluation of pi,j,i≥Y with ease for service time distributions with closed form Laplace transforms. The second expression of pi,j in Theorem 3 is similar to the case of the M/G/1 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 chainXn,n=0,⋯, the limiting probabilities πi,i=0,⋯,M−1 can be evaluated by solving the global balance equations πP=π;∑i=0M−1π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,
(27)πjpj,j−1=∑k=0j−1πk∑i=jM−1pk,i,which implies
(28)πj=∑k=0j−1πkak,jpj,j−1,where ak,j=∑i=jM−1pk,i.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(29)pi=limt⟶∞PXt=i,i=0,⋯,M.Moreover, we are interested in determining the performance measures of the system.LetES=1/μ be node 2 mean service time, λi=Nλ for i≤Y, and λi=M−iλ otherwise. Then, the throughput λf is given by
(30)λf=Nλπ0+NλES.The system size probabilities,pi,i=0,⋯,M, are given by
(31)pi=λfπiλi=πiπ0+NλES,i=0,⋯Y,NπiM−iπ0+NλES,i=Y+1,⋯M−1,and pM=1−∑i=0M−1pi.Moreover, the mean number of node 2 customersL is given by L=∑i=1Mipi. Little’s law implies that node 2 waiting time is W=L/λf. Furthermore, the mean delay is Wq=W−ES, 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.


## 4.1. Departure-Time Probabilities

Having determined the transition probability matrix for the Markov chainXn,n=0,⋯, the limiting probabilities πi,i=0,⋯,M−1 can be evaluated by solving the global balance equations πP=π;∑i=0M−1π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,
(27)πjpj,j−1=∑k=0j−1πk∑i=jM−1pk,i,which implies
(28)πj=∑k=0j−1πkak,jpj,j−1,where ak,j=∑i=jM−1pk,i.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(29)pi=limt⟶∞PXt=i,i=0,⋯,M.Moreover, we are interested in determining the performance measures of the system.LetES=1/μ be node 2 mean service time, λi=Nλ for i≤Y, and λi=M−iλ otherwise. Then, the throughput λf is given by
(30)λf=Nλπ0+NλES.The system size probabilities,pi,i=0,⋯,M, are given by
(31)pi=λfπiλi=πiπ0+NλES,i=0,⋯Y,NπiM−iπ0+NλES,i=Y+1,⋯M−1,and pM=1−∑i=0M−1pi.Moreover, the mean number of node 2 customersL is given by L=∑i=1Mipi. Little’s law implies that node 2 waiting time is W=L/λf. Furthermore, the mean delay is Wq=W−ES, 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.


## 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 from0 to infinity. Specifically, we select the deterministic, Erlang, exponential, and the hyperexponential distribution functions. Let bt be the p.d.f. of the service times. Recall that we require the service time distribution function to have a mean EB=1/μ. Here, we give explicit forms for the derivatives of the LST of the service time distribution functions.Deterministic. In this case, we assume that bt=b w.p. 1, so that
(32)Bn∗s=bne−sb.Exponential. Here bt=μe−μt,t≥0, so that
(33)Bn∗s=n!μs+μn+1.Erlang. The density function for bt=μμtk−1/k−1!e−μt,μ>0,t≥0. Here,
(34)Bn∗s=n!k+n−1k−1μks+μk+n.Hyperexponential. The two-phase hyperexponential is a mixture of two exponential distribution functions and has a density function of the form bt=2p2μe−2pμt+21−p2μe−21−pμt,μ>0,t≥0,0≤p≤1/2. Therefore,
(35)Bn∗s=n!2pμs+2pμn+1+21−pμs+21−pμn+1.For the above distributions, we recursively compute the stationary distributionsπ. and p. for the examples in the next subsection.


### 5.1. Small-Scale Applications

Here, we present a small-scale, two-node network with a total ofM=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.
(i)
For Theorem3 (i), where i=1,⋯,Y−1;j=i−1,⋯,Y, (i.e., i=1,2 and j=i−1,⋯,3), we use
(36)pi,j=bj−i+1=3λi−j+1Bi−j+1∗3λi−j+1!.(ii)
For Theorem3 (ii), where i=Y,⋯,M−1, j=i−1,⋯,M−1, (i.e., i=4,5 and j=2,3,4,5), we use
(37)pi,j=6−i!5−j!∑r=0j−i+1−1rB∗5−j+rλr!j−i+1−r!k!.(iii)
For Theorem3 (iii), where i=1,⋯,Y−1;j=Y+1,⋯,M−1(i.e., i=1,2 and j=4,5), we use the convolution result
(38)pi,j=34−i∑k=1j−2−1j−k−23−kCk,j−2a5−jk4−iB∗3−kλ−∑r=03−ikλrBr∗3λr!,where
(39)C1,1a=1,andC2,2a=2.Furthermore, fori=1,2 and j=5, we simply use pi,5=1−∑n=04pi,n. Moreover, pi,j=0 where j<i−1.Wheni=0, we use p0,j=p1,j for all j=0,⋯,M−1. Intuitively, if i is the number of customers at node 2 immediately after a departure, then the transition probabilities when i=0 must be the same as the probabilities when i=1, 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
(40)P=B∗3λ3λB1∗3λ4.5λ2B2∗3λ4.5λ3B2∗3λ47.5B∗3λ+6.75B∗λ−54B∗2λ+40.5λB1∗3λ+13.5λ2B2∗3λ1−∑j=04p0,jB∗3λ3λB1∗3λ4.5λ2B2∗3λ4.5λ2B3∗3λ47.5B∗3λ+6.75B∗λ−54B∗2λ+40.5λB1∗3λ+13.5λ2B2∗3λ1−∑j=04p1,j0B∗3λ3λB1∗3λ4.5λ2B2∗3λ4.5B∗λ−18B∗2λ+13.5B∗3λ+9λB1∗3λ1−∑j=04p2,j00B∗3λ3B∗2λ−3B∗3λ3B∗λ−6B∗2λ+3B∗3λ1−∑j=04p3,j000B∗2λ2B∗λ−2B∗2λ1−∑j=04p4,j0000B∗λ1−B∗λ.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,bt=b, and 0 otherwise (i.e., μ−1=b); thus, we have Bn∗s=bne−sb; therefore, B∗kλ=e−kλb; and Bn∗kλ=bne−kλb. This gives the following: B∗0=1,B∗λ=e−λb,B∗2λ=e−2λb,B∗3λ=e−3λb,B1∗3λ=be−3λb,B2∗3λ=b2e−3λb, and B3∗3λ=b3e−3λb. 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 tobw.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.
(41)P=e−3λb3λbe−3λb4.5λ2b2e−3λb4.5λ3b2e−3λb47.5e−3λb+6.75e−λb−54e−2λb+40.5λbe−3λb+13.5λ2b2e−3λb1−∑j=04p0,je−3λb3λbe−3λb4.5λ2b2e−3λb4.5λ2b3e−3λb47.5e−3λb+6.75e−λb−54e−2λb+40.5λbe−3λb+13.5λ2b2e−3λb1−∑j=04p1,j0e−3λb3λbe−3λb4.5λ2b2e−3λb4.5e−λb−18e−2λb+13.5e−3λb+9λbe−3λb1−∑j=04p2,j00e−3λb3e−2λb−3e−3λb3e−λb−6e−2λb+3e−3λb1−∑j=04p3,j000e−2λb2e−λb−2e−2λb1−∑j=04p4,j0000e−λb1−e−λb.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 ofM 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−PM−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 0,1,⋯,M−K+1. 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 (H2) distribution discussed earlier in this section. We fix the mean of the distribution to 1.0, i.e., μ=1, for all four distributions. For H2, 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.Table 1
The model withM=6, N=3, Y=3, and ρ=0.4.
D.F.p0p1p2p3p4p5p6λfLWDet.0.04530.10510.16840.24880.27330.13710.02190.95473.09873.2457Exp.0.10240.12290.14750.17700.21240.16990.06800.89763.05553.4042E20.07640.11920.15850.20430.23710.15960.04450.92363.06493.3185H20.26030.15900.09840.06220.04190.02550.35260.73972.95333.9927Table 2
The model withM=6, N=3, Y=3, and ρ=0.8.
D.F.p0p1p2p3p4p5p6λfLWDet.0.00010.00090.00760.06240.28760.46190.17960.99994.74054.7409Exp.0.00660.01580.03800.09110.21870.34990.27990.99344.66894.6998E20.00190.00710.02470.08450.24850.39600.23730.99814.70784.7166H20.05160.06280.07670.09390.11590.09890.50020.94844.45744.6999To 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 forM=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.Table 3
The model withM=100, N=80, Y=20, and ρ=0.012.
D.F.λf−AlgorithmL−AlgorithmL−SimulationDet..9959.1929.242E2.95010.91710.900Exp..94512.12612.113H2 (p=.005).74427.80827.884For the second example, we useM=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.Table 4
The model withM=60, N=40, Y=20, and ρ=0.025.
D.F.λf−AlgorithmL−AlgorithmL−SimulationDet..98013.06512.903E2.97213.61213.678Exp..96514.05814.165H2 (p=.125).93515.25615.317Note that changingp 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.Figure 3
Probability mass functions.Figure 4
Cumulative distribution functions.The hyperexponential presents the biggest challenge. To verify the correctness of our computations, we compared thecdf 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.


## 5.1. Small-Scale Applications

Here, we present a small-scale, two-node network with a total ofM=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.
(i)
For Theorem3 (i), where i=1,⋯,Y−1;j=i−1,⋯,Y, (i.e., i=1,2 and j=i−1,⋯,3), we use
(36)pi,j=bj−i+1=3λi−j+1Bi−j+1∗3λi−j+1!.(ii)
For Theorem3 (ii), where i=Y,⋯,M−1, j=i−1,⋯,M−1, (i.e., i=4,5 and j=2,3,4,5), we use
(37)pi,j=6−i!5−j!∑r=0j−i+1−1rB∗5−j+rλr!j−i+1−r!k!.(iii)
For Theorem3 (iii), where i=1,⋯,Y−1;j=Y+1,⋯,M−1(i.e., i=1,2 and j=4,5), we use the convolution result
(38)pi,j=34−i∑k=1j−2−1j−k−23−kCk,j−2a5−jk4−iB∗3−kλ−∑r=03−ikλrBr∗3λr!,where
(39)C1,1a=1,andC2,2a=2.Furthermore, fori=1,2 and j=5, we simply use pi,5=1−∑n=04pi,n. Moreover, pi,j=0 where j<i−1.Wheni=0, we use p0,j=p1,j for all j=0,⋯,M−1. Intuitively, if i is the number of customers at node 2 immediately after a departure, then the transition probabilities when i=0 must be the same as the probabilities when i=1, 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
(40)P=B∗3λ3λB1∗3λ4.5λ2B2∗3λ4.5λ3B2∗3λ47.5B∗3λ+6.75B∗λ−54B∗2λ+40.5λB1∗3λ+13.5λ2B2∗3λ1−∑j=04p0,jB∗3λ3λB1∗3λ4.5λ2B2∗3λ4.5λ2B3∗3λ47.5B∗3λ+6.75B∗λ−54B∗2λ+40.5λB1∗3λ+13.5λ2B2∗3λ1−∑j=04p1,j0B∗3λ3λB1∗3λ4.5λ2B2∗3λ4.5B∗λ−18B∗2λ+13.5B∗3λ+9λB1∗3λ1−∑j=04p2,j00B∗3λ3B∗2λ−3B∗3λ3B∗λ−6B∗2λ+3B∗3λ1−∑j=04p3,j000B∗2λ2B∗λ−2B∗2λ1−∑j=04p4,j0000B∗λ1−B∗λ.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,bt=b, and 0 otherwise (i.e., μ−1=b); thus, we have Bn∗s=bne−sb; therefore, B∗kλ=e−kλb; and Bn∗kλ=bne−kλb. This gives the following: B∗0=1,B∗λ=e−λb,B∗2λ=e−2λb,B∗3λ=e−3λb,B1∗3λ=be−3λb,B2∗3λ=b2e−3λb, and B3∗3λ=b3e−3λb. 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 tobw.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.
(41)P=e−3λb3λbe−3λb4.5λ2b2e−3λb4.5λ3b2e−3λb47.5e−3λb+6.75e−λb−54e−2λb+40.5λbe−3λb+13.5λ2b2e−3λb1−∑j=04p0,je−3λb3λbe−3λb4.5λ2b2e−3λb4.5λ2b3e−3λb47.5e−3λb+6.75e−λb−54e−2λb+40.5λbe−3λb+13.5λ2b2e−3λb1−∑j=04p1,j0e−3λb3λbe−3λb4.5λ2b2e−3λb4.5e−λb−18e−2λb+13.5e−3λb+9λbe−3λb1−∑j=04p2,j00e−3λb3e−2λb−3e−3λb3e−λb−6e−2λb+3e−3λb1−∑j=04p3,j000e−2λb2e−λb−2e−2λb1−∑j=04p4,j0000e−λb1−e−λb.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 ofM 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−PM−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 0,1,⋯,M−K+1. This reliability model is considered in example 4 by Maddah and El-Tha [1]. With our approach, we can solve larger reliability problems.


## 5.1.1. Deterministic Service Times

Let node 2 service times be a deterministic probability density function,bt=b, and 0 otherwise (i.e., μ−1=b); thus, we have Bn∗s=bne−sb; therefore, B∗kλ=e−kλb; and Bn∗kλ=bne−kλb. This gives the following: B∗0=1,B∗λ=e−λb,B∗2λ=e−2λb,B∗3λ=e−3λb,B1∗3λ=be−3λb,B2∗3λ=b2e−3λb, and B3∗3λ=b3e−3λb. 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 tobw.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.
(41)P=e−3λb3λbe−3λb4.5λ2b2e−3λb4.5λ3b2e−3λb47.5e−3λb+6.75e−λb−54e−2λb+40.5λbe−3λb+13.5λ2b2e−3λb1−∑j=04p0,je−3λb3λbe−3λb4.5λ2b2e−3λb4.5λ2b3e−3λb47.5e−3λb+6.75e−λb−54e−2λb+40.5λbe−3λb+13.5λ2b2e−3λb1−∑j=04p1,j0e−3λb3λbe−3λb4.5λ2b2e−3λb4.5e−λb−18e−2λb+13.5e−3λb+9λbe−3λb1−∑j=04p2,j00e−3λb3e−2λb−3e−3λb3e−λb−6e−2λb+3e−3λb1−∑j=04p3,j000e−2λb2e−λb−2e−2λb1−∑j=04p4,j0000e−λb1−e−λb.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 ofM 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−PM−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 0,1,⋯,M−K+1. 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 (H2) distribution discussed earlier in this section. We fix the mean of the distribution to 1.0, i.e., μ=1, for all four distributions. For H2, 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.Table 1
The model withM=6, N=3, Y=3, and ρ=0.4.
D.F.p0p1p2p3p4p5p6λfLWDet.0.04530.10510.16840.24880.27330.13710.02190.95473.09873.2457Exp.0.10240.12290.14750.17700.21240.16990.06800.89763.05553.4042E20.07640.11920.15850.20430.23710.15960.04450.92363.06493.3185H20.26030.15900.09840.06220.04190.02550.35260.73972.95333.9927Table 2
The model withM=6, N=3, Y=3, and ρ=0.8.
D.F.p0p1p2p3p4p5p6λfLWDet.0.00010.00090.00760.06240.28760.46190.17960.99994.74054.7409Exp.0.00660.01580.03800.09110.21870.34990.27990.99344.66894.6998E20.00190.00710.02470.08450.24850.39600.23730.99814.70784.7166H20.05160.06280.07670.09390.11590.09890.50020.94844.45744.6999To 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 forM=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.Table 3
The model withM=100, N=80, Y=20, and ρ=0.012.
D.F.λf−AlgorithmL−AlgorithmL−SimulationDet..9959.1929.242E2.95010.91710.900Exp..94512.12612.113H2 (p=.005).74427.80827.884For the second example, we useM=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.Table 4
The model withM=60, N=40, Y=20, and ρ=0.025.
D.F.λf−AlgorithmL−AlgorithmL−SimulationDet..98013.06512.903E2.97213.61213.678Exp..96514.05814.165H2 (p=.125).93515.25615.317Note that changingp 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.Figure 3
Probability mass functions.Figure 4
Cumulative distribution functions.The hyperexponential presents the biggest challenge. To verify the correctness of our computations, we compared thecdf 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.



---

*Source: 5550363-2021-05-30.xml*