Distribution function of the shortest path in networks of queues

OR Spectrum (2005) 27: 123–144 DOI: 10.1007/s00291-004-0169-3 c Springer-Verlag 2005
Distribution function of the shortest path in networks of queues Amir Azaron1 and Mohammad Modarres2 1 2 Department of Industrial Engineering, University of Bu-Ali-Sina, Ghobare Hamadani Boulevard, Hamadan, Iran (e-mail: aazaron@msl.sys.hiroshima-u.ac.jp) Department of Industrial Engineering, Sharif University of Technology, Azadi Street, Tehran, Iran (e-mail: modarres@sharif.edu) Abstract. In this paper, we consider acyclic networks of queues as a model to support the design of a dynamic production system. Each service station in the network represents a manufacturing or assembly operation. Only one type of product is produced by the system, but there exist several distinct production processes for manufacturing this product, each one corresponding with a directed path in the network of queues. In each network node, the number of servers in the corresponding service station is either one or infinity. The service time in each station is either exponentially distributed or belongs to a special class of Coxian distribution. Only in the source node, the service system may be modeled by an M/G/∞ queue. The transport times between every pair of service stations are independent random variables with exponential distributions. In method proposed in this paper, the network of queues is transformed into an equivalent stochastic network. Next, we develop a method for approximating the distribution function of the length of the shortest path of the transformed stochastic network, from the source to the sink node. Hence, the method leads to determining the distribution function of the time required to complete a product in this system (called the manufacturing lead time). This is done through solving a system of linear differential equations with non-constant coefficients, which is obtained from a related continuous-time Markov process. The results are verified by simulation. Keywords: Queueing – Stochastic processes – Graph theory – Shortest path Correspondence to: M. Modarres 124 A. Azaron and M. Modarres 1 Introduction Networks of queues are one of the most important issues in queueing theory, due to their vast application in a variety of domains. The often complex new problems that arise from these application domains are a source of inspiration for many researchers. A network of queues contains nodes that represent queueing systems or stations, each consisting of either one or an infinite number of parallel servers. The pattern of service for each customer can be unique. In the field of production or service management, many problems are modeled in the frame of a network of queues. An example is a dynamic production system, which produces only one type of product. An acyclic network of queues can represent this production system, in which the service stations are manufacturing departments. The total time each product is in a department is equal to the processing plus waiting time. If we assume several production processes are possible in principle, then each path in the network corresponds with one type of production process. Demand arrives at the service station settled in the source node according to a Poisson process while each finished product leaves the system from the sink node. The transport times between every two service stations are also modeled as independent random variables with exponential distributions. The length of each path in the network is equal to the sum of the lengths of the arcs and the nodes of this path. The length of each node is equal to the sojourn time in the particular queueing system (waiting time in queue plus the processing time), and the arc lengths are equal to the transport times between the service stations. Hence, the manufacturing lead time (or the total time spent by a product in the production system) is equivalent to the length of the shortest path in the network of queues. In this paper, we develop a method to obtain the manufacturing lead time distribution of a dynamic production system, which, in our model, turns into the development of an approach to approximate the steady-state distribution function of the shortest path in a network of queues. In this regard, Azaron and Fatemi Ghomi [1] developed a model for optimal control of service and arrival rates of a special class of Jackson networks, in order to minimize the expected shortest length of the network, as well as the average operating cost of the service stations. In each service station, the number of servers is either one or infinite, while the service time is a random variable with either an exponential distribution or a distribution from a special class of Coxian distribution functions (under heavy traffic conditions). Only the service time of the queue located in the source node may have any general distribution, while the number of servers in that node is assumed to be infinite (hence, an M/G/∞ system, since arrivals follow a Poisson process). In practice, a queueing system with infinite servers indicates that there is ample capacity so that no product ever has to wait. In the method proposed in this paper, we first transform each network of queues into a stochastic network. Then, the distribution function of the shortest path of this stochastic network (from the source node to the sink node) is determined through solving a system of differential equations with non-constant coefficients. By applying a technique from continuous-time Markov processes, this system of Distribution function of the shortest path in networks of queues 125 differential equations is constructed. To verify the method, the results are compared against the simulation results. Shortest path analysis in networks of queues has not been studied extensively so far; to our knowledge, no paper can be found that focuses on shortest path calculations. However, several methods are proposed in the literature for finding the shortest path in stochastic or dynamic networks. Several authors have addressed the problem, analytically. Martin [17] developes a method to obtain the distribution function of shortest path in stochastic networks as well as its expected value, where the arc lengths are independent random variables with polynomial distributions in the form of multiple integrals. Frank [7] computes the probability that the duration of the shortest path in a stochastic network is less than a specific value, when the arc lengths are continuous random variables. Kulkarni [14] presents an exact method to compute the distribution function of the length of the shortest path in directed stochastic networks, in which the arc lengths are independent random variables with exponential distributions. He constructs a continuous-time Markov chain with a single absorbing state and proves that the time until absorption into this state is equal to the length of the shortest path in the original network provided it starts from the initial state. Corea and Kulkarni [5] propose a methodology for computing the distribution function of the length of the shortest path, as well as criticality indices of paths. They assume that the arc lengths are integer valued and construct Markov chains with absorbing states and binary transition costs. Hayhurst and Shier [11] propose a method that is based on structural factoring. Sigal et al. [21] use uniformly directed cuts in their analysis of shortest path in directed acyclic networks. The intractability of these problems in general has motivated some researchers to develop approximation methods. Mirchandani [18] derives a bound on the expected length of a shortest path. By assuming that the arc lengths are discrete random variables, he avoids the complexity of numerical evaluation of multiple integrals. Hassin and Zemel [10] derive a bound for the expected length of a shortest path for uniformly distributed arcs with a fixed mean. Lyons et al. [16] discuss a lower bound for the expected shortest path length, based upon so-called resistance between the nodes on that path, or in other words, the expected length of the arcs connecting these nodes. Preiditis and Davis [19] derive an exact summation formula and a closed form approximation for the expected length of the shortest path for a complete graph (a graph which has an arc between every two nodes), where the arc lengths are independent and exponentially distributed random variables, based on Kulkarni’s method [14]. Kleindorfer’s method [12] also relies on bounds. Due to the difficulties arising in exact computation schemes, some researches have focused on Monte Carlo simulations for estimating the distribution function of the length of a shortest path (Sigal et al. [22]; Fishman [6]). Burt and Garman [4] develop a simulation method by conditional Monte Carlo simulation and the use of antithetical random variables for proper variance reduction techniques. They obtain good estimations of the distribution of the length of the shortest path in stochastic networks. A few publications have addressed the problem of determining shortest path in stochastic dynamic networks. Hall [9] introduces the problem of finding the least 126 A. Azaron and M. Modarres expected travel time path between two nodes in a network, in which travel times are random and time-dependent. Psaraftis and Tsitsiklis [20] examine shortest path problems, in which the arc costs are known functions of certain environmental variables at network nodes, and each of these variables evolves according to an independent Markov process. Vehicles traveling through the network can wait at a node in anticipation of more favorable arc costs. Azaron and Kianfar [2] generalize the paper by Psaraftis and Tsitsiklis [20]. They assume the environmental variables to evolve according to some independent semi-Markov process, rather than a Markov process. They also assume the length of each arc to be an exponential random variable, of which the parameter is a function of the environmental variable of its initiative node. Reviewing the above papers indicates that no work can be found similar to the current research devoted to determining a shortest path analysis in networks of queues. Although we develop an approximation method to obtain the distribution function of shortest paths, no bounds or Monte Carlo simulation approaches are applied. When the network does not contain any M/Ck /1 queueing system, our method is even exact. In fact, we generalize the method developed by Kulkarni [14] in order to find the distribution of the manufacturing lead time in a special class of dynamic production systems. This paper is organized in the following way. Section 2 reviews the shortest path analysis in stochastic networks, briefly. In Section 3, shortest path analysis in networks of queues is presented. We develop a shortest path algorithm in Section 4. In Section 5, the method is illustrated by some numerical examples as well as by a comparison of the results against those of simulation experiments. Finally, we draw conclusions in Section 6. 2 Shortest path analysis in stochastic networks In this section, we review Kulkarni’s method [14] to obtain the distribution function of the shortest path length (from the source to the sink node) of a directed acyclic stochastic network. We assume arc lengths are mutually independent random variables with exponential distributions. Let G = (V, A) be a directed acyclic stochastic network, where V and A represent the sets of nodes and arcs of the network, respectively. Let s and y represent the source and the sink node of the network, respectively. Let l (u, v) be the length of arc (u, v) ∈ A, which is an exponential random variable with parameter λ(u,v) . To construct a proper stochastic process, it is convenient to visualize the stochastic network as a communication network, in which a node is considered to be a station capable of receiving and transmitting messages and arcs as one-way communication links connecting pairs of nodes. The messages are assumed to travel at unit speed so that l (u, v) denotes the travel time from node u to v. As soon as a node receives a message over an incoming arc, it transmits it along all outgoing arcs and then disables itself (i.e., loses the ability to receive and transmit any future messages). This process continues until the message reaches the sink node y. Now, at any time there may be some nodes and arcs in the stochastic network that are “useless” for the progress of the message towards the sink node, i.e., even if the Distribution function of the shortest path in networks of queues 127 2 1 4 5 3 Fig. 1. The stochastic network corresponding with Example 1 messages are received and transmitted by these nodes and carried by these arcs, the message can only reach disabled nodes. It is assumed that all such “useless” nodes are also disabled and that messages traveling on such arcs are aborted. Now, let X(t) be the set of all disabled nodes at time t. X(t) is called the state of the network at time t. Definition 1. For each subset of nodes X ⊂ V , where s ∈ X and y ∈ X = V −X, we define the following sets: 1. X 1 = {v ∈ X: each path that connects any node of this set to the sink node y, must visit at least one member of X}. 2. S(X) = X ∪ X 1 . Definition 2. The state space Ω ∗ of the stochastic process to be constructed is defined by Ω = {X ⊂ V |X = S(X)}, (1) Ω ∗ = Ω ∪ {V }. (2) Definition 3. If X ⊂ V such that s ∈ X and y ∈ X, then a cut is defined as:
 C X, X = {( u, v) ∈ A/u ∈ X, v ∈ X}. (3) 
There is a unique minimal cut contained in C X, X , denoted by C(X). If X ∈ Ω, then
 C X, X = C(X). Example 1. Consider the network of Figure 1. Let X = (1, 2), then X 1 = φ, and S(X) = (1, 2). However, if X = (1, 3, 4), then the only path that connects node (2)∈ X to node (5) visits node (4), which belongs to X. Therefore, X 1 = (2), and S(X) = (1, 2, 3, 4). In this case, node (2) and arcs (1,2) and (2,4) would be “useless”. If X = (1, 2, 4), then node (3) does not belong to X 1 because it can be connected to (5) directly and the path 3-5 does not visit any node of X. In this case, X 1 = φ, and S(X) = (1, 2, 4). In this example, Ω ∗ = {(1), (1, 2), (1, 3), (1, 2, 3), (1, 2, 4), (1, 2, 3, 4), (1, 2, 3, 4, 5)}. The first six elements of Ω ∗ are the members of Ω while the last element of this set is V . If l (u, v) s, (u, v) ∈ A, are independent exponentially distributed random variables, then {X(t), t ≥ 0} is a continuous-time Markov process with state space Ω ∗ 128 A. Azaron and M. Modarres Table 1. The generator matrix Q for the example network State 1 1 2 3 4 5 6 7 0 0 0 0 0 0 2 3 4 5 6 7 λ(1,2) λ(1,3) 0 0 λ(1,3) λ(1,2) 0 λ(2,4) 0 0 0 0 λ(3,4) λ(2,4) + λ(3,4) λ(1,3) 0 0 λ(3,5) λ(3,5) λ(4,5) λ(3,5) + λ(4,5) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 and the infinitesimal generator matrix Q = [q(X, Y )](X, Y ∈ Ω ∗ ), see Kulkarni [14] for details, where   λ(u,v) if Y = S(X ∪ {v}),      (u,v)∈C(X)  (4) [q(X, Y )] = − λ(u,v) if Y = X,    (u,v)∈C(X)   0 otherwise. Whenever there is a transition in the process {X(t), t ≥ 0}, the number of disabled nodes in the stochastic network increases by at least one. Let us now order the elements of Ω ∗ by non-decreasing cardinality, i.e., if D, B ∈ Ω ∗ and |D| < |B|, then D appears before B in this ordering. From now on we assume that the states in Ω ∗ are numbered 1, 2, . . ., N = |Ω ∗ | so that the generator matrix Q is obviously upper triangular. Thus state 1 represents the source node, while the absorbing state N denotes the set of all nodes V . In Example 1, state 1 is {1}, and state 7 is {1, 2, 3, 4, 5}. Table 1 shows the infinitesimal generator matrix Q for this example. Each diagonal element of the matrix Q is equal to the minus sum of the other elements in the particular row. Let T represent the length of the shortest path. Clearly, T = min{t > 0 : X(t) = N/X(0) = 1}. Thus the length of the shortest path in the stochastic network is equal to the time until {X(t), t ≥ 0} gets absorbed in the final state N , starting from state 1. The objective is to compute F (t) = P {T ≤ t} or the distribution function of the length of the shortest path in the stochastic network. The Chapman-Kolmogorov backward or forward equations can be applied to compute F (t). Applying the backward algorithm, if we define, Pi (t) = P {X(t) = N/X(0) = i}1 ≤ i ≤ N, (5) then clearly we have F (t) = P1 (t). Using the backward algorithm, the system of differential equations for the vector P (t) = [P1 (t), P2 (t), . . . , PN (t)]T is given by Distribution function of the shortest path in networks of queues 129 P ′ (t) = Q.P (t), P (0) = [0, 0, . . ., 1]T . (6) By taking advantage of the upper triangular nature of Q, we present an efficient algorithm to compute F (t) in the transformed stochastic network, in Section 4. 3 Shortest path analysis in networks of queues In this section, first we propose a method for transforming an acyclic network of queues (called the original network) into an equivalent stochastic network, by replacing each node containing a queueing system with a stochastic arc. The length of this arc is equal to the sojourn time in the queueing system located in the corresponding node of the original network. Next, we derive/approximate the steadystate sojourn time distributions of the queueing systems settled in the nodes of the original network, except for the source node, by exponentials. Finally, we compute the steady-state distribution function of the length of the shortest path in the transformed stochastic network, using a method to be outlined in Section 4. Let us explain how to replace node k of the network of queues, which contains a queueing system, with a stochastic arc. Assume that b1 , b2 , . . ., bn are the incoming arcs to this node and d1 , d2 , . . ., dm are the outgoing arcs from it. Then, we substitute this node by arc (k ′ , k ′′ ), whose length is equal to the sojourn time in system for the particular queueing system. Furthermore, all arcs bi for i = 1, . . ., n end up with k ′ while all arcs dj for j = 1, . . ., m start from node k ′′ . Therefore, the network of queues is transformed into a stochastic network. The indicated process is the opposite of the absorption of an edge e in a graph G in graph theory (G.e), see Bondy and Murty [3] for more details. Since the length of the arc in the transformed stochastic network is equal to the sojourn time of the queueing system located in the node of the original network, we need to derive the sojourn time distributions of different classes of queuing systems in the next sub-sections. 3.1 Networks of queues with M/M/∞ queueing systems In this case, the sojourn time in system is equal to the service time with exponential distribution, because there is no queue. Therefore, we can replace each node, which includes an M/M/∞ queueing system, with an exponential arc with a parameter equal to its service rate. 3.2 Networks of queues with M/M/1 queueing systems It is well-known (and easily derived) that the density function of the sojourn time of an M/M/1 queueing system is exponentially distributed with parameter (µ − λ), hence w(t) = (µ − λ)e−(µ−λ)t t > 0, (7) 130 A. Azaron and M. Modarres where λ and µ are the arrival rate and the service rate of this queueing system, respectively. Hence, each single server node is replaced by an exponential arc with parameter µ − λ). 3.3 Networks of queues with an M/G/∞ queueing system As mentioned earlier, only the source node may be an M/G/∞ system. The service distribution of this system may be any arbitrary continuous function with hazard rate function λ(t). Let f (t) and F (t) be the density function and the distribution function of service time, respectively. Then, λ(t) is given by λ(t) = f (t)/(1 − F (t)). (8) In this case, we should replace λ with λ(t) in the infinitesimal generator matrix (4) and consequently, this matrix is a function of t. After transforming the network of queues into a stochastic one, we can directly solve the system of linear differential equations with non-constant coefficients (9), using the forward algorithm, in order to compute the distribution function of the length of the shortest path, by P̂ ′ (t) = P̂ (t).Q(t), P̂ (0) = [1, 0, . . ., 0]T , (9) where F (t) = P̂N (t) (in Sect. 4, we present a more efficient algorithm to obtain the distribution function of shortest path in the transformed stochastic network). Although the density function of the service time in an M/G/∞ queueing system is not exponential, the steady-state density function of time between two successive departures from this queueing system is exponential, see e.g. Gross and Harris [8]. Therefore, the departure process from an M/G/∞ queueing system is a Poisson process. The rate of this process is equal to the arrival rate to this queueing system. We point out once more that the above equations for the M/G/∞ system are valid only in the source node of the network. Any other node having a general service distribution would have a hazard rate function equal to λ(x), where x is the amount of time that the product has already been in that station. This hazard rate function is not a function of t, and therefore the system of differential equations (9) is not applicable to compute the distribution function of the length of the shortest path. 3.4 Networks of queues with a special class of M/Ck /1 queueing systems Approximating general distributions by phase-type (P H) distributions has significant applications in the analysis of queueing systems. A popular approach to analyzing queueing systems involving a general distribution G is to approximate G by a P H distribution. A P H distribution is a very general mixture of exponential distributions. It has been shown that matching three moments is sufficient for the accurate modeling of almost any system. Most existing algorithms for fitting Distribution function of the shortest path in networks of queues 131 a general distribution G to a P H distribution are restricted to a subset of P H distributions, being the class of Coxian distributions. Now, consider a k-stage Coxian distribution, Ck . The residence time in stage i is exponentially distributed with mean 1/µi , and the Markov chain is entered with probability pi in stage i, i = 1, 2, . . ., k. In this paper, we restrict our attention to a special class of Coxian distributions, where pi = 1, i = 1, 2, . . ., k. This class of Coxian distributions is also known as the class of generalized Erlang distributions (in Erlang distributions, µ1 = µ2 = ... = µk = µ). Even in this case, computing the density function of sojourn time is not easy, compared to the previous cases. Therefore, we consider the situation where each M/Ck /1 queueing system satisfies heavy traffic conditions, meaning that its utilization factor approaches 1. Now, we compute the diffusion approximation of the density function of the waiting time in queue, which has a high accuracy in this special case, see e.g. Gross and Harris [8] or Kleinrock [13]. Theorem 1. If an M/Ck /1 queueing system, where pi = 1, i = 1, 2, . . ., k, satisfies heavy traffic conditions, then it can be replaced by k + 1 exponential serial arcs with parameters (−2γ)/σ 2 , µ1 , µ2 , ..., µk , where γ and σ 2 denote the expected value and the variance of the waiting time in queue, respectively. Proof. Considering heavy traffic assumption, the expected value and the variance of the waiting time in queue in each M/G/1 queueing system are γ∆t = (λE(S) − 1)∆t + o(∆t), σ 2 ∆t = λE(S 2 )∆t + o(∆t), (10) where S and λ represent the service time and the arrival rate, respectively. Then, the utilization factor is equal to ρ = λE(S). Let wq (x, t/x0 ) represent the density function of the waiting time in queue at time t, given that the queueing system has started from x0 . Then, wq (x, t/x0 ) can be computed by solving the Fokker-Planck equation ∂wq (x, t/x0 )/∂t = −γ∂q (x, t/x0 )/∂x + σ 2 /2∂ 2 wq (x, t/x0 )/∂x2 , with the following boundary conditions:  ∞ wq (x, t/x0 )dx = ρ, (11) 0 wq (x, t/x0 ) ≥ 0. (12) Since the network of queues is in steady-state, the limiting density function exists, i.e., wq (x) = limt→∞ wq (x, t/x0 ) and can be computed as follows: σ 2 /2d2 wq (x)/dx2 − γwq (x)/dx = 0. (13) This is a second order differential equation with constant coefficients. By considering the following limiting boundary conditions under heavy traffic conditions,  ∞ wq (x)dx = ρ ∼ = 1, 0 wq (x) ≥ 0, (14) 132 A. Azaron and M. Modarres we may conclude that 2 wq (x) = (−2γ)/σ 2 e(2γ)/σ x x > 0. (15) Therefore, the steady-state density function of the waiting time in queue is approximated by an exponential density function with parameter (−2γ)/σ 2 . We can decompose each random variable with this special class of Coxian distribution into a collection of structured exponentially distributed random variables with parameters µi , i = 1, 2, . . ., k. Since the distribution of the waiting time in queue is exponential with parameter (−2γ)/σ 2 and the service time can be decomposed into k exponential serial arcs with parameters µi , i = 1, 2, . . ., k, Theorem 1 is proved. ⊓ ⊔ 3.5 Networks of queues with G/M/1 queueing systems If the arrival to each queueing system follows a Poisson process, all queueing systems are Markovian. In that case, let rij represent the probability of movement of a product from node i to node j and let λi be the rate of arrival Poisson process to the service station settled in node i. Then,  rji λj i ∈ V. (16) λi = j∈V However, even if there is only one M/G/1 queueing system in the network, the subsequent queueing systems connected to this one no longer face Poisson arrival processes. We now assume that the departure process from an M/G/1 queueing system settled in node u ∈ V , has independent increments. This assumption is reasonable, because the network is acyclic. Therefore, we can consider the departure process from node u as a non-homogeneous Poisson process with an intensity function equal to the hazard rate function of the time between two successive departures from this queueing system. Let B(t) and ρ represent the distribution function of service time and the utilization factor of the M/G/1 queueing systems settled in every node u. Then, C(t) or the distribution function of time between two successive departures from this queueing system is given by  t −λu C(t) = ρB(t) + (1 − ρ) B(t − u)λe du. (17) 0 The rate of the departure process from the queueing system settled in node u or λu (t) is equal to the hazard rate function of the time between two successive departures from this queueing system. Therefore, λu (t) is given by λu (t) = C ′ (t)/(1 − C(t)). Consequently, N (t) or the departure process from the queueing system settled in node u would be a non-homogeneous Poisson process with intensity function λu (t), and with the following distribution:  t m t − λu (s)ds 0 m! m ≥ 0. (18) λu (s)ds P [N (t) = m] =e 0 Distribution function of the shortest path in networks of queues 133 Theorem 2. Let N1 (t) represent the arrival process from the queueing system settled in node u to the queueing system settled in node v ∈ V . Then, N1 (t) is a non-homogeneous Poisson process with its intensity function equal to ruv λu (t), hence with distribution  t n t − r λ (s)ds n! n ≥ 0. (19) ruv λu (s)ds P [N1 (t) = n] = e 0 uv u 0 Proof. We compute the distribution of N1 (t) by conditioning on N (t): ∞  P [N1 (t) = n] = P [N1 (t) = n/N (t) = m].P [N (t) = m]. (20) m=n P [N1 (t) = n/N (t) = m] is the probability of n successes and m − n failures in m independent trials, where ruv is the probability of success on each trial. That is,  m n P [N1 (t) = n/N (t) = m] = (1 − ruv )m−n . (21) ruv n N (t) is a non-homogeneous Poisson process with the distribution (18). Consequently, P [N1 (t) = n]  ∞  m − n = ruv (1 − ruv )m−n e n m=n t 0 λu (s)ds  m t λu (s)ds m! , 0 P [N1 (t) = n] t = n − ruv e ∞  0 λu (s)ds n t  n! . λu (s)ds 0  (1 − ruv )m−n m−n t (m − n)!. λu (s)ds 0 m=n ∞ m−n m=n (1 − ruv ) t (1−ruv ) e 0 λu (s)ds  t 0 λu (s)ds m−n  (m − n)! is obviously equal to , and consequently, P [N1 (t) = n] t − n = ruv e 0 λu (s)ds  t  t n (1−ruv ) n! .e λu (s)ds t 0 λu (s)ds , 0 P [N1 (t) = n] t − =e 0 ruv λu (s)ds 0 n ruv λu (s)ds n! n ≥ 0. ⊓ ⊔ Corollary 1. The rate of arrival process to each G/M/1 queueing system, settled in node k ∈ V , is  rjk λj (t) k ∈ V. (22) λk (t) = j∈V 134 A. Azaron and M. Modarres Proof. The interarrival time of the queueing system settled in node k is the minimum of the interarrival time of the queueing systems settled in every node j to the queueing system settled in node k, j ∈ V. The hazard rate function of this random variable is rjk λj (t). Therefore, the hazard rate function of the minimum of these random variables, which indicates  the rate of arrival process to the queueing system ⊓ ⊔ settled in node k, is equal to j∈V rjk λj (t). The closed form (22) is quite similar to (16). After computing λk (t), we can compute A(t), or the distribution function of the time between two successive arrivals to node k, from the following equation: t − A(t) = 1 − e 0 λk (u)du . (23) Let µ be the service rate of the G/M/1 queueing system, then we can compute 0 < x0 < 1, the unique root of equation (24), by using a numerical method like the Newton-Raphson method,  ∞ e−µt(1−z) dA(t) 0 < z < 1. (24) z= 0 After computing x0 , we can compute the density function of sojourn time in the G/M/1 queueing system, by: w(t) = µ(1 − x0 )e−µ(1−x0 )t t > 0. (25) Clearly, w(t) has exponential distribution with parameter µ(1 − x0 ), and we can replace the queueing system settled in this node by an exponential arc with parameter µ(1 − x0 ). 4 The shortest path algorithm In this section, we first present an efficient method to compute the distribution function of shortest path in the transformed stochastic network. If an M/G/∞ queueing system settles in the source node of the network of queues, then the first and second elements of the first row of the infinitesimal generator matrix Q(t) will have non-constant values, while all other elements are equal to zero. Therefore, we can directly solve the system of linear differential equations with non-constant coefficients (9) to obtain the distribution function of shortest path. Due to the upper triangular nature of Q(t) and the fact that there is at most one M/G/∞ queueing system settled in the source node of the network, we can develop a more efficient method to compute the distribution function of the length of the shortest path. Let Q represent an N − 1 × N − 1 sub-matrix of Q(t) with the elimination of the first row and the first column of Q(t). Now, let us explain how the system of differential equations with constant coefficients (26), is solved. If P (t) = [P2 (t), P3 (t), ..., PN (t)]T , then, ′ P (t) = Q.P (t), T P (0) = [0, 0, . . ., 1] . (26) Distribution function of the shortest path in networks of queues 135 Let M be the modal matrix of Q, an N − 1 × N − 1 matrix whose N − 1 columns are the eigenvectors of Q. Let α1 , α2 , ..., αN −1 be the eigenvalues of Q, which are the diagonal elements of Q owing to the upper triangular nature of Q. We can compute P (t) from equation (27), P (t) = M eAt M −1 P (0), (27) At where e is the diagonal matrix, see Luenberger [15] for the details, as follows:  α1 t  e 0 .0  0 eα2 t . . . eAt =   . . .. αN −1 t 0 . .e According to the definition of P2 (t), the distribution function of the length of the shortest path in the transformed stochastic network would be the convolution of P2 (t) and f (t) or the density function of service time in the M/G/∞ queueing system settled in the source node of the network of queues. Finally, the distribution function F (t) of the length of the shortest path in the network of queues is computed from the following equation  t P2 (t − s)f (s)ds. (28) F (t) = 0 If there is no M/G/∞ queueing system in the network of queues, then we can easily compute P (t) in the system of differential equations (6). Considering M as the modal matrix of Q, A as the diagonal matrix whose diagonal elements are the eigenvalues of Q and P (0) = [0, 0, . . ., 1]T , P (t) is computed from equation (29). P (t) = M eAt M −1 P (0). (29) In this case, F (t) = P1 (t). The complete algorithm for the computation of F (t) now reads as follows. Algorithm Step 1.1. Compute the arrival process rate to the M/M/1 and M/Ck /1 queueing systems from equation (16). Step 1.2. Compute the arrival process rate to each G/M/1 queueing system, settled in node k ∈ V , from equation (22). Step 2.1. Compute the steady-state distribution of the sojourn time in system for each M/M/1 queueing system from equation (7) with λ as determined in Step 1.1. Step 2.2. For the M/M/∞ and M/G/∞ queueing systems, the sojourn time distribution is equal to the service time distribution. Step 2.3.1. For each G/M/1 queueing system, compute A(t) from equation (23) with λk (t) as computed in Step 1.2. Step 2.3.2. Compute x0 , or the unique root of equation (24), numerically. Step 2.3.3. Compute the sojourn time distribution for each G/M/1 queueing system from equation (25). Step 3.1. Compute γ and σ 2 from equations (10). 136 A. Azaron and M. Modarres Step 3.2. Compute the distribution of the waiting time in queue, for each M/Ck /1 queueing system, from equation (15). Step 4. Replace each node containing an M/M/∞, M/G/∞, M/M/1 or G/M/1 queueing system with a stochastic arc whose length is equal to the sojourn time for the particular queueing system. Step 5. Replace each node containing an M/Ck /1 queueing system with k + 1 exponential serial arcs. The length of the first arc is equal to the waiting time in queue for this queueing system, taken from Step 3.2, and the lengths of the next serial arcs are exponentially distributed with the parameters µ1 , µ2 , ..., µk (parameters of the Coxian distribution of the service time). Step 6. Eliminate all arcs with zero length. Step 7. If the source node contains an M/G/∞ queueing system, then, go to Step 8.1. Otherwise, go to Step 9.1. Step 8.1. Compute M, the modal matrix of Q, and eAt (A is a diagonal matrix whose elements are the eigenvalues of Q). Step 8.2. Compute P2 (t) from equation (27). Step 8.3. Compute F (t) from equation (28). Then, go to Step 10. Step 9.1. Compute M , the modal matrix of Q, and eAt (A is a diagonal matrix whose elements are the eigenvalues of Q). Step 9.2. Compute F (t) = P1 (t) from equation (29). Step 10. Stop. The steady-state distribution function of shortest path in the network of queues, or the manufacturing lead time distribution in the dynamic production system is equal to F (t). 5 Numerical examples To illustrate the proposed approach and the results, we solve 5 numerical examples. Example I is depicted in Figure 2. This dynamic production system produces only one type of product. There are two production processes to produce this product, each one corresponding with a directed path of the network of queues. The product demand arrives at the queueing system settled in the source node according to a Poisson process with the rate of λ = 1 per hour and each finished product leaves the system from the sink node. The objective is to obtain the distribution function of the manufacturing lead time in this production system. Node A contains an M/G/∞ queueing system with a Weibull distribution 2 of the service time, with parameters (α, β) = (1, 2), f (t) = 2te−t t > 0. Node B contains an M/C2 /1 queueing system with the service parameters (µ1 , µ2 , p1 , p2 ) = (1, 2, 1, 1). Node C contains a G/M/1 queueing system with the service rate equal to µC = 4. Node D contains an M/M/1 queueing system with the service rate equal to µD = 2. There is no queueing system in node E. The lengths of arcs k, l, m and p are zero, but the length of arc n has exponential distribution with parameter 3. It is also assumed that rAB = 0.66 and rAD = 0.34. Distribution function of the shortest path in networks of queues 137 m B k A C D l p E n Fig. 2. The network of queues corresponding with the Examples I, III and V 2t 1 0.0086 2 1 3 1.66 2 5 3.7 6 7 3 4 Fig. 3. The stochastic network corresponding with Example I The arrival rate at the M/M/1 queueing system settled in node D is equal to 0.34. The arrival rate at the M/C2 /1 queueing system settled in node B is equal to 0.66 and its utilization factor is ρB = 0.99. The rate of arrival process to the G/M/1 queueing system settled in node C is equal to the rate of departure process from the M/C2 /1 queueing system settled in node B, which is also equal to the hazard rate function of the service time in this queueing system, because ρB approaches one. Therefore, A(t), i.e., the distribution function of time between two successive arrivals to the queueing system settled in node C is equal to the service time distribution of M/C2 /1 queueing system. Then, we compute x0 , or the unique root of the following equation.  ∞ z= e−4t(1−z) (2e−t − 2e−2t )dt0 < z < 1, 0 x0 = 0.075. (30) The sojourn time distributions for the M/M/1 and G/M/1 queueing systems are exponential with the parameters equal to 1.66 and 3.7, respectively. We can also easily compute γ = −0.01 and σ 2 = 2.31. The distribution of the waiting time in queue for the M/C2 /1 queueing system is exponential with the parameter equal to 0.0086. The hazard rate function of the sojourn time of the M/G/∞ queueing system settled in node A is equal to 2t. Next, we replace the nodes A, C and D with the proper stochastic arcs. We also replace node B containing the M/C2 /1 queueing system with 3 exponential serial arcs with parameters equal to 0.0086, 1 and 2, respectively. Then, we eliminate the arcs k, l, m and p from the network. Finally, the network of queues is transformed into a proper stochastic network, depicted in Figure 3. Now, we find the distribution function of shortest path from node 1 to node 7 in this stochastic network (the number over each arc shows the hazard rate function of the length of particular arc). The indicated stochastic process {X(t), t ≥ 0} has 10 states in the following order: Ω ∗ = {(1), (1, 2), (1, 2, 3), (1, 2, 4), (1, 2, 3, 4), (1, 2, 3, 5), (1, 2, 3, 4, 5), (1, 2, 3, 5, 6), (1, 2, 3, 4, 5, 6), (1, 2, 3, 4, 5, 6, 7)}. Table 2 shows the infinitesimal generator matrix Q(t). 138 A. Azaron and M. Modarres Table 2. Matrix Q(t) corresponding with Example I 1 State 1 2 3 4 5 6 7 8 9 10 2 3 4 5 6 7 8 9 10 −2t 2t 0 0 0 0 0 0 0 0 0 −1.6686 0.0086 1.66 0 0 0 0 0 0 0 0 −2.66 0 1.66 1 0 0 0 0 0 0 0 −3.0086 0.0086 0 0 0 0 3 0 0 0 0 −4 0 1 0 0 3 0 0 0 0 0 −3.66 1.66 2 0 0 0 0 0 0 0 0 −5 0 2 3 0 0 0 0 0 0 0 −5.36 1.66 3.7 0 0 0 0 0 0 0 0 −6.7 6.7 0 0 0 0 0 0 0 0 0 0 According to Steps 8.1, 8.2 and 8.3 of the proposed algorithm, F (t) results as follows:  t 2 2 2 2se−s + 0.0044se−s +5.36s−5.36t − 0.0024se−s +6.7s−6.7t F (t) = 0 2 2 − 0.042se−s +3.66s−3.66t + 0.0232se−s +5s−5t 2 2 + 0.1066se−s +2.66s−2.66t − 0.059se−s +4s−4t − 4.58se−s 2 +1.6686s−1.6686t + 2.5344se−s 2 +3.0086s−3.0086t  ds . (31) We can easily compute the moments of the distribution, using a suitable numerical integration technique. We use MATLAB 5.3 to compute the expected value and the variance of the manufacturing lead time. Theses values are as follows:  ∞ µ= tdF (t) = 1.8277, 0 σ2 =  ∞ t2 dF (t) − µ2 = 0.6839. (32) 0 We also solve 4 other numerical examples. All of these networks has 4 service stations, but with different configurations and characteristics of the service stations. All arc lengths are equal to zero. Example II is depicted in Figure 4. This is a serial network of queues. The demand arrives at the source node according to a Poisson process with the rate of λ = 1 per hour. Node A contains an M/G/∞ queueing system. The distribution of service time in this system is Weibull with parameters (α, β) = (1, 2). Nodes B and D contain the M/M/1 queueing systems with service rates of µB = µD = 2. Node C contains an M/M/∞ queueing system with the service rate equal to µC = 1. The hazard rate function of the sojourn time of the M/G/∞ queueing system settled in node A is equal to 2t. The exponential parameters of the sojourn times of the M/M/1, M/M/∞ and M/M/1 queueing systems settled in nodes B, C and D, respectively, are all equal to 1. Distribution function of the shortest path in networks of queues A B C 139 D Fig. 4. The network of queues corresponding with Example II 2t 1 2 1 1 1 3 5 4 Fig. 5. The transformed stochastic network corresponding with Example II 3 1 1 2 3 3 4 1 Fig. 6. The transformed stochastic network corresponding with Example III Figure 5 shows the transformed stochastic network. The number above each arc shows the hazard rate function of the length of particular arc. In this example, the size of the state space is equal to 5, Ω ∗ = {(1), (1, 2), (1, 2, 3), (1, 2, 3, 4), (1, 2, 3, 4, 5)}, and F (t) is given by  t 2 2 2 F (t) = 2se−s − s(t − s)2 e−s +s−t − 2s(t − s)e−s +s−t 0 − 2se−s 2 +s−t ds. (33) The configuration of the network of queues in Example III is the same as Example I, Figure 2. This network of queues does not have any non-Markovian queueing system. The arrival process is a Poisson process with the rate of λ = 2 per hour. Node A contains an M/M/∞ queueing system with the service rate equal to µA = 3. Nodes B, C and D contain the M/M/1 queueing systems with the service rates equal to µB = 2, µC = 4 and µD = 2, respectively. There is no queueing system in node E. It is also assumed that rAB = rAD = 0.5. The exponential parameters of the sojourn times of the M/M/∞, M/M/1, M/M/1 and M/M/1 queueing systems settled in nodes A, B, C and D, respectively, are equal to 3, 1, 3 and 1, respectively. Figure 6 shows the transformed stochastic network. The number over each arc shows the hazard rate function of its length. In this example, the size of the state space is equal to 4, Ω ∗ = {(1), (1, 2), (1, 2, 3), (1, 2, 3, 4)}, and F (t) is given by F (t) = 1 − 1.5e−4t − 4.5e−2t + 5e−3t . (34) Example IV is depicted in Figure 7, which has more than one G/M/1 queueing system. The arrival rate of Poisson process is λ = 0.49 per hour. Node A contains an M/G/∞ queueing system with Weibull distribution of service time and parameters (α, β) = (2, 2). Node B contains an M/C2 /1 queueing system with service parameters (µ1 , µ2 , p1 , p2 ) = (1, 1, 1, 1) (clearly, this is an M/E2 /1 queueing system with service parameters (µ, k) = (1, 2)). Nodes C and D contain G/M/1 queueing systems with service rates equal to µC = 4 and µD = 2, respectively. 140 A. Azaron and M. Modarres A B C E D Fig. 7. The network of queues corresponding with Example IV 4t 1 0.014 2 1 3 1 4 3.885 5 1.956 6 Fig. 8. The transformed stochastic network corresponding with Example IV There is no queueing system in node E. It is also assumed that rBC = 0.8 and rBD = 0.2. The hazard rate function of the sojourn time of the M/G/∞ queueing system settled in node A is 4t. The unique root of equation (24) for the G/M/1 queueing systems settled in nodes C and D are equal to 0.0287 and 0.0221, respectively. Therefore, the exponential parameters of the sojourn of the G/M/1 queueing systems settled in nodes C and D will be equal to 3.885 and 1.956, respectively. The exponential parameter of the waiting time in queue for the M/E2 /1 queueing system settled in node B is equal to 0.014. The exponential parameters of the next serial arcs corresponding with the service time of this queueing system are both equal to 1. Figure 8 shows the transformed stochastic network. The number over each arc shows the hazard rate function of its length. In this example, the size of the state space is equal to 6, Ω ∗ = {(1), (1, 2), (1, 2, 3), (1, 2, 3, 4), (1, 2, 3, 4, 5), (1, 2, 3, 4, 5, 6)}, and F (t) is given by  t 2 2 F (t) = 3.88se−2s − 4se−2s +0.014s−0.014t 0 + 0.068s(t − s)e−2s 2 +s−t + 0.12se−2s 2 +s−t ds. (39) The configuration of the network of queues in Example V is the same as Example I, Figure 2. This is the first network of queues that we analyze and has more than one M/Ek /1 queueing system. The arrival process to the queueing system settled in the source node of this network is assumed to be a Poisson process with the rate of λ = 0.98 per hour. Node A contains an M/G/∞ queueing system with the Weibull distribution of service time and the parameters (α, β) = (1, 2). Nodes B and D contain both the M/E2 /1 queueing systems with service parameters (µ, k) = (1, 2). Node C contains a G/M/1 queueing system with the service rate equal to µC = 4. There is no queueing system in node E. It is also assumed that rAB = rAD = 0.5. The hazard rate function of the sojourn time of the M/G/∞ queueing system settled in node A is equal to 2t. The value of x0 for the G/M/1 queueing system settled in node C is equal to 0.0429. Therefore, the exponential parameter of the sojourn time of the G/M/1 queueing system settled in node C would be equal to 3.828. The exponential parameters of the waiting time in queue for the M/E2 /1 Distribution function of the shortest path in networks of queues 2t 0.014 2 1 0.014 1 141 1 3 3.828 8 6 5 1 1 4 7 Fig. 9. The transformed stochastic network corresponding with Example V 1.2 1 0.8 ,H VD0.6 & F(t) Fsim(t) 0.4 0.2 0       W      Fig. 10. F (t) and F sim(t) versus t, in Example I queueing systems settled in nodes B and D are both equal to 0.014. The exponential parameters of the next serial arcs corresponding with the service times of these queueing systems are all equal to 1. Figure 9 shows the transformed stochastic network. The number over each arc shows the hazard rate function of the length of particular arc. In this example, the size of the state space is equal to 14, Ω ∗ = {(1), (1, 2), (1, 2, 3), (1, 2, 4), (1, 2, 3, 4), (1, 2, 3, 5), (1, 2, 4, 7), (1, 2, 3, 4, 7), (1, 2, 3, 5, 6), (1, 2, 3, 4, 5), (1, 2, 3, 4, 5, 7), (1, 2, 3, 4, 5, 6), (1, 2, 3, 4, 5, 6, 7), (1, 2, 3, 4, 5, 6, 7, 8)}, and F (t) is given by  t 2 2 2 F (t) = 1.998se−s − 0.002se−s +0.014s−0.014t + 0.004se−s +s−t 0 + 0.002s(t − s)e−s 2 +s−t − 2se−s + 0.001se−s 2 +3.842s−3.842t − 0.002se−s 2 +2s−2t 2 +0.028s−0.028t + 0.066s(t − s)e−s − 0.002s(t − s)e−s  2 − 0.0006s(t − s)2 e−s +2s−2t ds . 2 2 +1.014s−1.014t +2s−2t (40) For verifying the results, we use the Monte Carlo simulation to estimate the distribution functions of the shortest path in the networks of queues of Examples I, IV and V, F sim(t), (the values of F (t) in Examples II and III are exact). The sample size is n = 1000. Figures 10, 11 and 12 show F (t) and F sim(t) versus time for Examples I, IV and V, respectively. According to these figures, the maximum difference between F (t) and F sim(t), in these three examples, is equal to 2.97%. Therefore, comparing the simulation results against our approximation results shows the accuracy of our proposed method. 142 A. Azaron and M. Modarres 1.2 1 0.8 ,9 HV0.6 D& F(t) Fsim(t) 0.4 0.2 0                  W Fig. 11. F (t) and F sim(t) versus t, in Example IV 1.2 1 0.8 9 H VD0.6 & F(t) Fsim(t) 0.4 0.2 0                    W Fig. 12. F (t) and F sim(t) versus t, in Example V 6 Conclusion In this paper, we have developed an approximation method for computing the steady-state distribution function of shortest path in the networks of queues. This method is constructed based on the theory of stochastic processes, queueing theory and graph theory. Furthermore, we assumed the arc lengths to be independent random variables with exponential distributions. The importance of this subject arises from its vast application in various real world problems such as in production systems, reliability modeling and computer networks. Many problems are formulated in the framework of network of queues and consequently this method can be applied. For example, we obtained the distribution function of the manufacturing lead time in a class of dynamic production systems. In the proposed method, the network of queues was transformed into a proper stochastic network. Then, we presented an efficient method for computing the distribution function of the length of the shortest path in the transformed stochastic network. Validity of the results was verified by simulation. In the case of general service time distribution, it is possible to approximate the sojourn time distribution by an appropriate Coxian distribution. To do that, we Distribution function of the shortest path in networks of queues 143 match the first three moments and then apply our proposed method to compute the distribution function of shortest path in the network of queues. For further research, it is necessary to understand the limitations of our method and then to develop some approaches to overcome the problems, which result from the limitations. The major problem facing this method arises from its nature of NP hardness. The state space of the continuous-time Markov process grows exponentially with the network size. From the examples presented in Section 5, it is obvious that the size of the state space depends on the number of service stations, number of arcs as well as the number of M/Ck /1 queueing systems. As the worst case example, for a completely transformed network with n nodes and n(n − 1) arcs, the size of the state space would be 2n−2 + 1. One must also note that for very large networks every method of producing reasonably accurate answers will be prohibitively expensive. In this paper, we analyzed the acyclic networks with no direct or indirect feedback or overtaken. Therefore, considering the sojourn times as independent random variables is a reasonable assumption. Another suggestion for extension of this research is to relax the independence assumption of sojourn times of queueing systems settled in the nodes of the network. In reality, this assumption is not always true. This model can also be extended for the cases containing other types of queueing systems or for computing the distribution function of longest path in the networks of queues. References 1. Azaron A, Fatemi Ghomi SMT (2003) Optimal control of service rates and arrivals in jackson networks. European Journal of Operational Research 147: 17–31 2. Azaron A, Kianfar F (2003) Dynamic shortest path in stochastic dynamic networks: ship routing problem. European Journal of Operational Research 144: 138–156 3. Bondy J, Murty U (1976) Graph theory with applications. Elsevier, Amsterdam 4. Burt J, Garman M (1971) Conditional Monte Carlo: a simulation technique for stochastic network analysis. Management Science 18: 207–217 5. Corea GA, Kulkarni VG (1993) Shortest paths in stochastic networks with ARC lengths having discrete distributions. Networks 23: 175–183. 6. Fishman GS (1985) Estimating network characteristics in stochastic activity networks. Management Science 31: 579–593 7. Frank H (1969) Shortest paths in probabilistic graphs. Operations Research 17: 583–599 8. Gross D, Harris M (1985) Fundamentals of queueing theory, 2nd Wiley, New York 9. Hall R (1986) The fastest path through a network with random time-dependent travel times. Transportation Science 20: 182–188 10. Hassin R, Zemel E (1985) On shortest paths in graphs with random weights. Mathematics of Operations Research 10: 557–564 11. Hayhurst KJ, Shier DR (1991) A factoring approach for the stochastic shortest path problem. Operations Research Letters 10: 329–334 12. Kleindorfer GB (1971) Bounding distributions for a stochastic acyclic network. Operations Research 19: 1586–1601 13. Kleinrock L (1975) Queueing systems, 2 Computer applications. Wiley, New York 14. Kulkarni VG (1986) Shortest paths in networks with exponentially distributed arc lengths. Networks 16: 255–274 144 A. Azaron and M. Modarres 15. Luenberger D (1979) Introduction to dynamic systems. Wiley, New York 16. Lyons R, Pemantle R, Peres Y (1999) Resistance bounds for first-passage percolation and maximum flow. Journal of Combinatorial Theory Ser A 86: 158–168 17. Martin J (1965) Distribution of time through a directed acyclic network. Operations Research 13: 46–66 18. Mirchandani P (1976) Shortest distance and reliability of probabilistic networks. Computers and Operations Research 3: 347–355 19. Preiditis A, Davis R (1993) The expected length of a shortest path. Information Processing Letters 46: 135–141 20. Psaraftis H, Tsitsiklis J (1993) Dynamic shortest paths in acyclic networks with markovian arc costs. Operations Research 41: 91–101 21. Sigal C, Pritsker A, Solberg J (1979) The use of cut sets in monte carlo analysis of stochastic networks. Mathematics and Computer Simulation 21: 376–384 22. Sigal C, Pritsker A, Solberg J (1980) The stochastic shortest route problem. Operations Research 28: 1122–1129