Analysis of Queueing Networks with Blocking using a New Aggregation Technique

Johann Christoph Strelen, Berthold Bärk, Jürgen Becker, Volker Jonas Rheinische Friedrich–Wilhelms–Universität Bonn

Abstract

Queueing networks with repetitive-service blocking are investigated. The state probabilities are calculated applying a new aggregation technique, the disaggregation-aggregation (DA) iteration. Thus one obtains close approximate solutions with product form.

Introduction

Queueing networks have been successfully used in modeling a variety of computer systems, communication networks, and manufacturing systems. They can efficiently be analyzed if they have a product form solution. In queueing networks with finite capacities, blocking occurs ([11, 12]). Only in special cases they have a product form solution ([1]), otherwise the analysis is difficult due to complexity.

We investigate queueing network models with repetitive-service blocking: If a job tries to enter a node, and there is no buffer space available, the job goes back to its former node and instantaneously receives another service there. Afterwards, a new target node is chosen. In this paper, the steady state probabilities are approximated applying disaggregation-aggregation iteration ([13, 14]).

Disaggregation aggregation (DA) iteration is a new general approximate aggregation technique for finite Markov chains with discrete time. The calculated state probabilities are in product form. Here, neither the complete state space, nor all state probabilities, nor the transition probability matrix are handled explicitly, marginal probabilities are considered instead. Therefore Markov chains with really large state spaces can be treated, for example billions of states. It is also possible to calculate time-dependent state probabilities.

Time and space demands are higher than for many known, special purpose heuristics, but are by orders of magnitude lower than for exact solutions of Markov chains with many states. On the other hand, the results of DA iteration are usually more accurate than the heuristic solutions. Moreover, the accuracy can be influenced: The more aggregates are considered, the more accuracy can be obtained, if the aggregation is chosen suitably. If many aggregates are used, the aggregation is called fine.

This work was supported in part by the Deutsche Forschungsgemeinschaft

September 8, 1997

It is interesting to note that the steady state probabilities of separable queueing networks are equal to DA solutions calculated according to suitably chosen aggregates.

We analyzed quite a lot of Markov chain models applying DA iteration, and we compared the obtained results with the exact solutions and with results which were calculated using heuristic methods known from the literature. The accuracy of DA solutions compares favourably.

We discuss DA iteration in section 1, its application to blocking networks in section 2, and numerical examples in section 3.

1 Disaggregation-aggregation iteration

Irreducible aperiodic Markov chains (Z(t), t ∈ IN), discrete and homogeneous in time, with finite state space Z = [1 : n], n ∈ IN, are considered. They stem from stochastic models of computer, communication, or manufacturing systems. These systems consist of interacting subsystems, their components. At every moment, each component is in one of its possible component states. This is why we denote the whole system’s states also by state tuples z = (z_{1}, …, z_{K}), where z_{k} indicates the state of component k. The states of component k are numbered z_{k} = 0, …, N_{k} for all components k = 1, …, K.

Aggregates are sets of system states. An aggregate can be defined by the state of a certain component, for example the aggregate Z_{k,j} is the set of all system states where the component k is in state j, irrespective of the states of the other components. All aggregates Z_{k,0}, …, Z_{k,N_{k}} form a partition of the state space, for k fixed.

In the sequel, we will use the notions “partition” and “aggregate” slightly more generally. Therefore we introduce them more formally: A class {Z_{k,0}, …, Z_{k,N_{k}}}, N_{k} ∈ IN, k ∈ IN, of disjoint non-empty subsets of Z whose union is the full set Z itself is called partition Z_{k} of Z. A class {Z_{1}, …, Z_{C}}, C ∈ IN, of partitions is called an aggregation AGG of Z, and each of the subsets Z_{k,j} of system states is an aggregate.

State tuples z = (z_{1}, …, z_{C}) ∈ [0 : N_{1}] × … × [0 : N_{C}] indicate the aggregates Z_{k,z_{k}} of each partition Z_{1}, …, Z_{C} to which a system state belongs. If for all state tuples the intersections Z_{1,z_{1}} ∩ … ∩ Z_{C,z_{C}} contain at most one element, each state tuple defines at most one system state z ∈ Z, and the aggregation is said to completely identify the states. If the intersection is empty, there is no according state, and the state tuple z is termed infeasible, otherwise z describes a state and is feasible. In the sequel we consider only aggregations which completely identify the states. Thus, for each state z ∈ Z there is one and only one such state tuple z. The set Z^{T} of all these feasible state tuples is another representation of the state space Z. There is a one-to-one mapping I : Z^{T} onto Z.

A special aggregation is defined by the structure of the model as follows. Consider a system which consists of K subsystems, for example a queueing network with K nodes. If each aggregate is defined by the state z_{k} of a single subsystem, Z_{k,j} = {I(z) | z ∈ Z^{T} , z_{k} = j}, j = 0, …, N_{k}, k = 1, …, K, where [0 : N_{k}] is the set of all states of subsystem k, then we call the aggregation canonical. If the whole system’s state is uniquely defined by the states of all subsystems, this canonical aggregation completely identifies the states. Consider for example a closed queueing network with N jobs circulating in it. z_{k} denotes the number of jobs in node k. We assume that these numbers z_{k} describe the states of the nodes completely. Under these circumstances, any state tuple (z_{1}, …, z_{K}) where z_{1} + … + z_{K} = N is feasible and indicates one and only one system state. If z_{1} + … + z_{K} ≠ N, the tuple is infeasible – the state (z_{1}, …, z_{K}) cannot occur.

p^{(t)}(z) =^{△} P{Z(t) = z}, z ∈ Z, denotes the state probabilities of the Markov chain at time t. For convenience, we use also p^{(t)}(z), z ∈ Z^{T}, where p^{(t)}(z) = p^{(t)}(z) holds if z and z denote the same state, i.e. z = I(z) holds. In the sequel, if z and z occur in the same context, they denote the same state.

The probability vector p^{(t)} = [p^{(t)}(1), …, p^{(t)}(n)] contains all state probabilities at time t, and P =^{△} [p(x, y)]_{x,y=1}^{n} denotes the matrix of the transition probabilities. For convenience, we use also p(x, y), x, y ∈ Z^{T}, where p(x, y) = p(x, y) for x = I(x), y = I(y) holds.

Sometimes, if we are not interested in a particular time, or if steady state probabilities are considered, we omit the specification (t).

Each aggregate Z_{k,j} defines a macro probability or aggregate probability which is a marginal probability as follows:

r_{k}(j) =^{△} P{Z ∈ Z_{k,j}} = ∑_{z∈Z_{k,j}} p(z), j = 0, …, N_{k}, k = 1, …, C.  (1)

We collect the marginal probabilities into tuples r_{k} =^{△} (r_{k}(0), …, r_{k}(N_{k})), k = 1, …, C,

and, in turn, these into R =^{△} (r_{1}, …, r_{C}).

By (1), the aggregation function

A : [0, 1]^{n} → [0, 1]^{N_{1}+1} × … × [0, 1]^{N_{C}+1}, p ↦ R,

is defined. It calculates the macro probabilities for given state probabilities. Sometimes joint marginal probabilities

r_{k,k′}(j, j′) =^{△} P{Z ∈ Z_{k,j}, Z ∈ Z_{k′,j′}} = ∑_{z∈Z_{k,j}∩Z_{k′,j′}} p(z),

j = 0, …, N_{k}, j′ = 0, …, N_{k′}, {k, k′} ⊆ [1 : C],

are used. Again we collect them into tuples r_{k,k′}.

Each partition Z_{k}, k ∈ [1 : C], defines a random variable Z_{k} : Z^{T} → [0 : N_{k}], z ↦ z_{k}. The marginal distribution r_{k} is the probability distribution of Z_{k}. Furthermore, sometimes it is more convenient to consider the random vector Z = (Z_{1}, …, Z_{C}) for which P{Z = z} = p(z), z ∈ Z^{T}, holds, instead of the random variable Z.

Our proposed solution technique needs disaggregation: Given positive macro probabilities R, a distribution π = [π(1), …, π(n)] over Z is defined as follows. First, it satisfies the macro probabilities, R = A(π), and secondly, π maximizes the entropy. Thus, according to the product form theorem, see Kagan, Linnik, Rao [7, theorem 13.2.1, page 409], this maximum entropy distribution (MED) π is uniquely determined and has a product form,

π(z) = \frac{1}{G} ∏_{k=1}^{C} ρ_{k,z_{k}},  (2)

where the ρ_{k,z_{k}} > 0 are some real numbers,

r_{k}(j) = \frac{1}{G} ∑_{{z∈Z^{T} \atop z_{k}=j}} ∏_{i=1}^{C} ρ_{i,z_{i}}, j = 0, …, N_{k}, k = 1, …, C,  (3)

holds, and

G = ∑_{z∈Z^{T}} ∏_{k=1}^{C} ρ_{k,z_{k}}  (4)

is a normalizing constant. In fact, the MED is completely determined by all MED parameters ρ_{k,j}, ρ for short.

By (2), (3), and (4), the disaggregation function is defined:

D : [0, 1]^{N_{1}+1} × … × [0, 1]^{N_{C}+1} → [0, 1]^{n}, R ↦ π.

It calculates the MED over the state space for given macro probabilities.

In general, the nonlinear equations (3) and (4) must be solved in order to calculate all MED parameters ρ_{k,j} and the normalizing constant G of the MED π, given the marginal probabilities R. This is easily accomplished in some special cases. For example, if a canonical aggregation completely identifies the states, and if the state space is rectangular, i.e. if Z^{T} is a cartesian product, Z^{T} = [0 : N_{1}] × … × [0 : N_{C}], the MED is given by

p(z) = r_{1}(z_{1}) · . . . · r_{C}(z_{C}), z ∈ Z^{T} ,  (5)

which is easily verified (see [13]). Here the MED parameters are ρ_{k,j} = r_{k}(j), j = 0, …, N_{k}, k = 1, …, C, and G = 1 holds. This can be generalized for the case where Z^{T} consists of the union of cartesian products, [14, section 3]. Another simple case where the state space is not necessarily rectangular are the state probabilities in product form queueing networks.

Instead of calculating the steady state probabilities q of a Markov chain we determine an approximation p = D(R) where R = T (R) holds. Here the function T is defined as follows:

T : [0, 1]^{N_{1}+1} × … × [0, 1]^{N_{C}+1} → [0, 1]^{N_{1}+1} × … × [0, 1]^{N_{C}+1}, R ↦ A(D(R)P).

This approximation p and its marginal probabilities R are referred to as disaggregation aggregation product form (DAΠ) solution. In fact, only the marginal probabilities R and the MED parameters ρ and G must be calculated. This is practicable even if the state space is large, so large that neither all state probabilities nor the transition probabilities could be handled.

We calculate the DAΠ solutions using disaggregation aggregation (DA) iteration,

R^{(t+1)} = T (R^{(t)}), t = 1, 2, …  (6)

until max_{k,j} |r_{k}^{(t+1)}(j) − r_{k}^{(t)}(j)| is sufficiently small, beginning with feasible marginal probabilities R^{(1)}.

This DA iteration resembles vector iteration, p^{(t+1)} = p^{(t)}P, t = 1, 2, …, beginning with any distribution p^{(1)}. p^{(t)} are the time-dependent state probabilities and A(p^{(t)}) the time-dependent macro probabilities. p^{(t)} converges to the stationary probability vector and A(p^{(t)}) to the stationary macro probabilities. Accordingly, when DA iteration R^{(t+1)} = T (R^{(t)}), t = 1, 2, …, is performed, R^{(t)} are approximations of the timedependent macro probabilities and D(R^{(t)}) approximations of the time-dependent state probabilities. If DA iteration is done for a sufficient number of iterations, R^{(t)} is an approximation of the stationary macro probabilities, and D(R^{(t)}) an approximation of the stationary state probabilities.

It is important to note that the function T must be evaluated without handling either all transition probabilities or all probabilities π^{(t)}. Whether this is possible, depends on the transition probabilities and on the aggregation. According to our experience, such conditions are fulfilled very often in Markov chains originating from performance or reliability models. We have not yet found any model which could not be treated with DA iteration although we tried many different models. In [14, section 3], conditions are given under which this is possible.

More precisely, the conditions are fulfilled if the following expressions for the joint probabilities φ_{k}(j, l) of the random vectors (Z_{k}^{(t)}, Z_{k}^{(t+1)}),

φ_{k}(j, l) = ∑_{{y∈Z^{T} \atop y_{k}=j}} ∑_{{z∈Z^{T} \atop z_{k}=l}} p(y, z)p^{(t)}(y), j, l = 0, …, N_{k}, k = 1, …, C,  (7)

can be simplified substantially by collecting many summands into only a few expressions. Observing

r_{k}^{(t+1)}(l) = ∑_{j=0}^{N_{k}} φ_{k}(j, l), l = 0, …, N_{k},  (8)

φ_{k}(j, j) = r_{k}^{(t)}(j) − ∑_{{l=0 \atop l≠j}}^{N_{k}} φ_{k}(j, l), j = 0, …, N_{k}, k = 1, …, C,  (9)

follows. (9) is already simple in the sense mentioned above.

The conditional probabilities

π_{k}(j, l) =^{△} P\{Z_{k}^{(t+1)} = l | Z_{k}^{(t)} = j\} = \frac{φ_{k}(j, l)}{r_{k}^{(t)}(j)}, j, l = 0, …, N_{k}, k = 1, …, C, (10)

are aggregate transition probabilities as the following aggregate transition equations

r_{k}^{(t+1)}(l) = ∑_{j=0}^{N_{k}} π_{k}(j, l)r_{k}^{(t)}(j), l = 0, …, N_{k}, k = 1, …, C,  (11)

hold. The aggregate transition probabilities are not homogeneous in time, in general. According to (9), the aggregate non-transition probabilities are

π_{k}(j, j) = 1 − ∑_{{l=0 \atop l≠j}}^{N_{k}} π_{k}(j, l), j = 0, …, N_{k}, k = 1, …, C.  (12)

(8) and (9) or (11) and (12) are used to efficiently calculate the new macro probabilities R^{(t+1)} from the old ones, R^{(t)}, thus evaluating the function T .

If no closed form solution is available, the nonlinear equations (3) and (4) can be solved iteratively as follows ([6]) for the MED parameters ρ, given the macro probabilities R:

ρ_{k,j}^{[i+1]} = (1 − ω_{i} + ω_{i} \frac{r_{k}(j)}{r_{k}^{[i]}(j)})ρ_{k,j}^{[i]}, j = 0, …, N_{k}, k = 1, …, K, i = 1, 2, …,  (13)

where the damping factor is chosen as ω_{i} = max{1, 1/2, 1/4, …} so that ρ^{[i+1]} fulfils the equations (3) and (4) more accurate than ρ^{[i]}. Here ρ^{[1]} are initial approximations of the unknown parameters ρ, and r_{k}^{[i]}(j) are the macro probabilities according to the MED parameters ρ^{[i]}, see (3). But in fact, the iteration (13) can be omitted when steady state solutions are desired. It is sufficient to apply just one iteration step,

ρ_{k,j}^{(t+1)} = (1 − ω_{t} + ω_{t} \frac{r_{k}^{(t+1)}(j)}{r_{k}^{(t)}(j)})ρ_{k,j}^{(t)}  (14)

within each DA iteration step (6). This iteration did never diverge, using a damping factor ω_{i} = 1 or ω_{i} = 1/2.

We summarize, stating the

DA Iteration Algorithm

Input: Reasonable MED parameters ρ^{(1)} – for example all ρ_{k,j}^{(1)} = 1 – and a small bound ϵ > 0

0) Calculate the normalization constant G^{(1)} and the marginal probabilities R^{(1)}. Let t = 1

1) Calculate the joint distribution functions φ_{k}(j, l) for all j, l, k, or the aggregate transition probabilities π_{k}(j, l) for all j, l, k

2) From these, calculate the new macro probabilities R^{(t+1)}

3) Calculate the MED parameters ρ^{(t+1)} either directly, for example using (5), or approximately using (14). In the latter case, calculate again the macro probabilities R^{(t+1)} from the new MED parameters ρ^{(t+1)}. Calculate the normalizing constant G^{(t+1)}

4) If max_{j,k} |r_{k}^{(t+1)}(j) − r_{k}^{(t)}(j)| < ϵ holds, the iteration is completed. Otherwise let t := t + 1 and goto step 1)

Results: The macro probabilities R^{(t+1)}, the MED parameters ρ^{(t+1)}, and the normalizing constant G^{(t+1)}. □

In [14] the DA product form theorem is proven: The equation R = T (R) has a solution.

Unfortunately, we have not yet found error estimates for DA solutions or a proof of convergence of the iteration. But in [14] is proven that there is always an aggregation with only a small number of aggregates for which an approximating product form (2) exists which is as accurate as desired. And we solved many models using DA iteration, among others polling, fork join, stochastic Petri net, and blocking models, and we found very accurate results; the iteration did never diverge.

Moreover, the accuracy can be substantially ameliorated if more and suitably chosen aggregates are considered. We tried this with a lot of models. An aggregation with many aggregates is called a fine grain aggregation.

But how should the aggregates be chosen? Two observations are important for this question. If a system consists of subsystems which behave stochastically independent, and if a Markov model for this system is solved using DA iteration with canonical aggregation, canonical with respect to these components, the DA solution is exact. On the other hand, if a system consists of highly dependent components, the canonical approximation is expected to be inaccurate. Consequently, accurate approximations can be expected if the components are nearly independent. This observation serves as a guide for aggregation: Look for independent components, and design the aggregation accordingly, or combine the clusters of dependent components, each cluster into a single new conmponent, in order to get nearly independent components. In general, this produces more aggregates and requires more computing time. These matters will be treated in [2].

In section 2, fine grain aggregations for queueing networks with blocking are presented, and in section 3, numerical results are presented which verify that accurate approximations can be obtained using fine grain aggregation.

With DA iteration, time-dependent state probabilities can be calculated as follows: If the Markov chain is discrete in time, the DA solutions p^{(t)}, t = 1, 2, …, are approximate time-dependent state probabilities, and R^{(t)} the according time-dependent macro probabilities. If the Markov chain is continuous in time, randomization ([5]) can be applied onto the DA solutions in order to obtain time-dependent solutions.

2 Queueing networks with repetitive-service blocking

We consider queueing networks with repetitive service blocking. A network consists of K ∈ IN nodes. In node k, k = 1, …, K, there is a single exponential server, service rate μ_{k}. Buffer space is available for no more than ν_{k} jobs. If the model is open, jobs arrive according to a Poisson process with parameter λ_{k}. When a service of a job in node k is finished, the job tries to go to node l with transfer probability t_{k,l}, l ∈ [1 : K]. If there is no buffer space available, the job remains in node k and instantaneously receives another service (repetitive-service with random destination policy). In an open network, a job which is completely served in node k leaves the network with probability t_{k,0}. In a closed network, there are N ∈ IN jobs, and we assume without loss of generality ν_{k} ≤ N, k = 1, …, K.

The system’s state is given by the K−tuple z = (z_{1}, …, z_{K}) where z_{k} denotes the number of jobs in node k. In an open queueing network, the state space

Z^{T} = S_{1} × … × S_{K}

where S_{k} = [0 : ν_{k}], and in a closed queueing network

Z^{T} = {z ∈ S_{1} × … × S_{K} | z_{1} + … + z_{K} = N},  (15)

where

S_{k} = [τ_{k} : ν_{k}], τ_{k} = max (0, N − ∑_{j≠k} ν_{j}) :

A state τ_{k} > 0 in node k occurs if no buffer space is left in any other node. □

The probabilities B_{k} that a job is blocked in node k, the blocking probabilities, are given by

B_{k} = ∑_{{j=1 \atop j≠k}}^{K} t_{k,j}(P\{Z_{j} = ν_{j}\} − P\{Z_{k} = 0, Z_{j} = ν_{j}\}), k = 1, …, K.  (16)

2.1 Open queueing networks

Here we treat open queueing networks with repetitive-service blocking. First, we give the Markov chain model. Secondly, DA iteration is presented applying canonical aggregation. Thirdly, fine grain solutions with more aggregates are described; in this manner even more accurate approximations are obtained.

The transitions from state y ∈ Z^{T} into state z ∈ Z^{T} occur as follows (1_{k} denotes the K-tuple with all elements 0 but 1 in position k):

1. arrival of a job at node k from outside, z = y + 1_{k}, if y_{k} < ν_{k}, rate λ_{k}

2. service at node k is completed:

(a) the job leaves the queueing network, z = y − 1_{k}, if y_{k} > 0, rate μ_{k}t_{k,0}

(b) the job goes from node k to node l, l ≠ k, z = y−1_{k}+1_{l}, if y_{k} > 0 and y_{l} < ν_{l}, rate μ_{k}t_{k,l}.

Thus the sparse generator matrix Q = [q_{y,z}]_{y,z=1}^{n} is given by

q_{y,z} = \{ \begin{matrix}μ_{k}t_{k,0} & if y_{k} > 0, z = I(y − 1_{k}), \\ μ_{k}t_{k,l} & if y_{k} > 0, y_{l} < ν_{l}, z = I(y − 1_{k} + 1_{l}), \\ λ_{k} & if y_{k} < ν_{k}, z = I(y + 1_{k}), \\ − ∑n_{{x=1 \atop x≠y}} q_{y,x} & if y = z, \\ 0 & otherwise, \end{matrix}  (17)

y = I^{−1}(y), y, z = 1, …, n.

In the sequel, we consider the Markov chain with the transition probability matrix P = QΔ + I where Δ = 0.99/ max_{x} |q_{x,x}|. By the factor 0.99, numerical problems stemming from rounding errors are avoided. This Markov chain is discrete in time and has the same steady state probabilities.

Canonical aggregation

Here the number C of partitions equals the number K of nodes in the queueing network. Aggregate Z_{k,j} is the set of all states with j jobs in node k,

Z_{k,j} = {I(z) | z ∈ Z^{T} , z_{k} = j}, k = 1, …, K, j = 0, …, ν_{k}.  (18)

Disaggregation is effectuated according to (5). The aggregate transition probabilities are ^{∗}

π_{k}(j, l) = \{ \begin{matrix}α_{k} & if l = j − 1 ≥ 0, \\ 1 − β_{k} & if l = j = 0, \\ 1 − α_{k} − β_{k} & if 0 < l = j < ν_{k}, \\ 1 − α_{k} & if l = j = ν_{k}, \\ β_{k} & if l = j + 1 ≤ ν_{k}, \\ 0 & otherwise, \end{matrix} j, l = 0, …, ν_{k},  (19)

with

α_{k} = μ_{k}Δt_{k,0} + ∑_{{c=1 \atop c≠k}}^{K} μ_{k}Δt_{k,c}(1 − r_{c}(ν_{c})), β_{k} = λ_{k}Δ + ∑_{{c=1 \atop c≠k}}^{K} μ_{c}Δt_{c,k}(1 − r_{c}(0)),

k = 1, …, K.

Observing the product form probabilities (5), the blocking probabilities are given by

B_{k} = (1 − r_{k}(0))( ∑_{{j=1 \atop j≠k}}^{K} t_{k,j}r_{j}(ν_{j})), k = 1, …, K.

The number of aggregates is ∑_{k=1}^{K}(ν_{k} + 1).

Fine grain aggregation

Here we use more aggregates, and thus we may expect – and in fact we obtain – more accurate solutions. More specific, we define aggregates according to the state of two nodes. Thus the dependencies between two such nodes are completely taken into account. The partitions are denoted by Z_{(k,k′)}, k, k′ ∈ [1 : K], k ≠ k′, and the aggregates are

Z_{(k,k′),(l,l′)} = {I(z) | z ∈ Z^{T} , z_{k} = l, z_{k′} = l′}, l = 0, …, ν_{k}, l′ = 0, …, ν_{k′}.

Not all partitions Z_{(k,k′)}, k ≠ k′, are considered, but only for a subset K_{2} ⊂ [1 : K]^{2} \ {(1, 1), …, (K, K)} of node pairs.

Here, the macro probabilities are the joint marginal probabilities r_{k,k′}(j, j′). They are calculated corresponding to (11) by the aggregate transition equations

r_{k,k′}^{(t+1)}(l, l′) = ∑_{j=0}^{ν_{k}} ∑_{j′=0}^{ν_{k′}} π_{k,k′}((j, j′), (l, l′))r_{k,k′}^{(t)}(j, j′), (l, l′) ∈ S_{k} × S_{k′}, (k, k′) ∈ K_{2}. (20)

In general, the aggregate transition probabilities π_{k,k′}((j, j′), (l, l′)) differ for different choices of the set K_{2}.

Sometimes the dependencies between a singular node k^{∗} ∈ [1 : K] and the other nodes are of particular importance, for example in a star topology. Then we try center aggregation with the partitions

Z_{(k^{∗},k)}, k = 1, …, K, k ≠ k^{∗}, K_{2} = {(k^{∗}, k)|k = 1, …, K, k ≠ k^{∗}},

∗Proofs are omitted in this paper. They are straightforward and sometimes lengthy.

and the partition Z_{k^{∗}}. Node k^{∗} is called the center. The number of aggregates is (ν_{k^{∗}} +1) (1 + ∑K_{{c=1 \atop c≠k^{∗}}} (ν_{c} + 1)) , and the MED calculated from given marginal probabilities is

p(z) = r_{k^{∗}}(z_{k^{∗}}) ∏_{{c=1 \atop c≠k^{∗}}}^{K} \frac{r_{k^{∗},c}(z_{k^{∗}}, z_{c})}{r_{k^{∗}}(z_{k^{∗}})}, z ∈ Z^{T} .  (21)

The aggregate transition probabilities are

π_{k,k′}((j, j′), (l, l′)) =

μ_{k^{∗}}Δt_{k^{∗},0} + ∑_{{c=1 \atop c≠k^{∗},k′}}^{C} μ_{k^{∗}}Δt_{k^{∗},c}(1 − \frac{r_{k^{∗}c}(j, ν_{c})}{r_{k^{∗}}(j)}) if l = j − 1 ≥ 0, l′ = j′,

μ_{k′}Δt_{k′,0} + ∑_{{c=1 \atop c≠k^{∗},k′}}^{C} μ_{k′}Δt_{k′,c}(1 − \frac{r_{k^{∗}c}(j, ν_{c})}{r_{k^{∗}}(j)}) if l = j, l′ = j′ − 1 ≥ 0,

λ_{k^{∗}}Δ + ∑_{{c=1 \atop c≠k^{∗},k′}}^{C} μ_{c}Δt_{c,k^{∗}}(1 − \frac{r_{k^{∗}c}(j, 0)}{r_{k^{∗}}(j)})  if l = j + 1 ≤ ν_{k}, l′ = j′,

λ_{k′}Δ + ∑_{{c=1 \atop c≠k^{∗},k′}}^{C} μ_{c}Δt_{c,k′}(1 − \frac{r_{k^{∗}c}(j, 0)}{r_{k^{∗}}(j)})  if l = j, l′ = j′ + 1 ≤ ν_{k′},
μ_{k}Δt_{k,k′}  if l = j − 1 ≥ 0, l′ = j′ + 1 ≤ ν_{k′},
μ_{k′}Δt_{k′,k}  if l = j + 1 ≤ ν_{k}, l′ = j′ − 1 ≥ 0,
0  if |l + l′ − j − j′| > 1,

j, l = 0, …, ν_{k}, j′, l′ = 0, …, ν_{k′}, (j, j′) ≠ (l, l′), (k, k′) ∈ K_{2}.

Here and in the sequel, the aggregate non-transition probabilities π_{k,k′}((j, j′), (j, j′)) are calculated according to (12). Aggregate transition probabilities for the partition Z_{k^{∗}} are not needed, the marginal probabilities r_{k^{∗}}(j), j ∈ S_{k^{∗}}, are obtained by (22).

The blocking probabilities are obtained using (16) where the marginal probabilities are given by

r_{k}(j) = \{ \begin{matrix}∑_{j′=0}^{ν_{k∗}} r_{k^{∗}k}(j′, j) & if k ≠ k^{∗}, \\ ∑_{j′=0}^{ν_{k′}} r_{k^{∗}k′}(j, j′) & if k = k^{∗}, \end{matrix} k = 1, …, K, k′ ∈ [1 : K] \ \{k^{∗}\}, (22)

r_{k,k′}(j, j′) = \{ \begin{matrix}∑_{j′′=0}^{ν_{k∗}} \frac{r_{k^{∗},k}(j′′, j)r_{k^{∗},k′}(j′′, j′)}{r_{k^{∗}}(j′′)} & if k ≠ k^{∗}, k′ ≠ k^{∗}, \\ r_{k^{∗}k′}(j, j′) & if k = k^{∗}, k′ ≠ k^{∗}, \\ r_{k^{∗}k}(j′, j) & if k ≠ k^{∗}, k′ = k^{∗}, \end{matrix}

j = 0, …, ν_{k}, j′ = 0, …, ν_{k′}, k, k′ = 1, …, K, k ≠ k′.

Sometimes the states of separate pairs (k, k′) of nodes depend strongly. Under such circumstances one can put these pairs into K_{2} and consider the according partitions Z_{(k,k′)}, (k, k′) ∈ K_{2}. If K is even, and all nodes of the network are combined into such pairs,

K_{2} = {(per(1), per(2)), (per(3), per(4)), …, (per(K − 1), per(K))}

holds, where (per(1), …, per(K)) is a permutation of (1, …, K). We call this paired aggregation.

The number of aggregates is ∑_{(k,k′)∈K_{2}}(ν_{k} + 1)(ν_{k′} + 1), and the MED is

p(z) = ∏_{(k,k′)∈K_{2}} r_{k,k′}(z_{k}, z_{k′}), z ∈ Z^{T} .

Here the aggregate transition probabilities are

π_{k,k′}((j, j′), (l, l′)) =

μ_{k}Δt_{k,0} + ∑_{{c=1 \atop c≠k,k′}}^{C} μ_{k}Δt_{k,c}(1 − r_{c}(ν_{c})) if l = j − 1 ≥ 0, l′ = j′,

μ_{k′}Δt_{k′,0} + ∑_{{c=1 \atop c≠k,k′}}^{C} μ_{k′}Δt_{k′,c}(1 − r_{c}(ν_{c})) if l = j, l′ = j′ − 1 ≥ 0,

λ_{k}Δ + ∑_{{c=1 \atop c≠k,k′}}^{C} μ_{c}Δt_{c,k}(1 − r_{c}(0))  if l = j + 1 ≤ ν_{k}, l′ = j′,

λ_{k′}Δ + ∑_{{c=1 \atop c≠k,k′}}^{C} μ_{c}Δt_{c,k′}(1 − r_{c}(0))  if l = j, l′ = j′ + 1 ≤ ν_{k′},
μ_{k}Δt_{k,k′}  if l = j − 1 ≥ 0, l′ = j′ + 1 ≤ ν_{k′},
μ_{k′}Δt_{k′,k}  if l = j + 1 ≤ ν_{k}, l′ = j′ − 1 ≥ 0,
0  if |l + l′ − j − j′| > 1,

j, l = 0, …, ν_{k}, j′, l′ = 0, …, ν_{k′}, (j, j′) ≠ (l, l′), (k, k′) ∈ K_{2}.

The blocking probabilities are obtained using (16) where the required marginal probabilities are given by

r_{k}(j) = \{ \begin{matrix}∑_{j′′=0}^{ν_{k}} r_{k′′,k}(j′′, j) & if (k′′, k) ∈ K_{2}, \\ ∑_{j′′=0}^{ν_{k′′}} r_{k,k′′}(j, j′′) & if (k, k′′) ∈ K_{2}, \end{matrix} j = 0, …, ν_{k}, k = 1, …, K,

r_{k,k′}(j, j′) = \{ \begin{matrix}r_{k,k′}(j, j′) & if (k, k′) ∈ K_{2}, \\ r_{k}(j)r_{k′}(j′) & otherwise, \end{matrix}
j = 0, …, ν_{k}, j′ = 0, …, ν_{k′}, k, k′ = 1, …, K.

Clearly, the states of two nodes k and k′ are dependent if the nodes are connected, i.e. if jobs travel from node k to node k′ or vice versa. Therefore it seems reasonable to consider all these pairs:

K_{2} = {(k, k′)|k < k′, t_{k,k′} > 0 or t_{k′,k} > 0}.

We call this the all connected pairs (ACP) aggregation.

First, we restrict ourselves to tandem queues where jobs arrive from outside at node 1 with rate λ_{1}, travel from node k to node k + 1, k = 1, …, K − 1, and leave the network at node K. Thus, K_{2} = {(1, 2), (2, 3), …, (K − 1, K)} is the set of considered pairs for ACP aggregation.

There are C = K −1 partitions of the state space Z^{T} , and the aggregate probabilities are joint marginal probabilities, namely

r_{k}(j, j′) =^{△} P{Z_{k} = j, Z_{k+1} = j′}, j ∈ S_{k}, j′ ∈ S_{k+1}, k = 1, …, C.

The number of aggregates is ∑_{c=1}^{C}(ν_{c} + 1)(ν_{c+1} + 1). Disaggregation is obtained by

p(z) = \frac{r_{1}(z_{1}, z_{2}) · . . . · r_{C}(z_{C}, z_{C+1})}{r_{2}(z_{2}) · . . . · r_{C}(z_{C})}, z ∈ Z^{T} .

From this, one obtains

r_{k}(j) = ∑_{j′=0}^{ν_{k+1}} r_{k}(j, j′), k = 1, …, C, and r_{k}(j) = ∑_{j′=0}^{ν_{k−1}} r_{k−1}(j′, j), k = 2, …, K.

The aggregate transition probabilities are

π_{k,k+1}((j, j′), (l, l′)) =

λ_{1}Δ  if k = 1, j + 1 = l, j′ = l′, j < ν_{1},
μ_{k−1}Δ(1 − r_{k−1}(0, j)/r_{k}(j))  if 1 < k, j + 1 = l, j′ = l′, j < ν_{k},
μ_{k}Δ  if j − 1 = l, j′ + 1 = l′, j > 0, j′ < ν_{k+1},
μ_{k+1}Δ(1 − r_{k+1}(j′, ν_{k+2})/r_{k+1}(j′)) if k < C, j = l, j′ − 1 = l′, j′ > 0,
μ_{K}Δ  if k = C, j = l, j′ − 1 = l′, j′ > 0,
0  otherwise if (j, j′) ≠ (l, l′),

j, l = 0, …, ν_{k}, j′, l′ = 0, …, ν_{k+1}, k = 1, …, C.

The blocking probabilities are obtained using (16), B_{k} = r_{k+1}(ν_{k+1}) − r_{k}(0, ν_{k+1}), k = 1, …, K − 1.

Now we discuss an aggregation which is even finer than the all connected pairs of nodes aggregation, the all connected triplets (ACT) aggregation. Here we restrict ourselves to open tandem queues as above. The triplets of nodes are (k, k +1, k +2), k = 1, …, K −2. There are C = K − 2 partitions Z_{k} of the state space with the aggregates

Z_{k,(j,j′,j′′)} = {I(z) | z ∈ Z^{T} , z_{k} = j, z_{k+1} = j′, z_{k+2} = j′′},
(j, j′, j′′) ∈ S_{k} × S_{k+1} × S_{k+2}, k = 1, …, C.

We will use the marginal probabilities r_{k}(j), the joint marginal probabilities r_{k}(j, j′), and the aggregate probabilities are the type-3 marginal probabilities

r_{k}(j, j′, j′′) =^{△} P{Z_{k} = j, Z_{k+1} = j′, Z_{k+2} = j′′}, (j, j′, j′′) ∈ S_{k}×S_{k+1}×S_{k+2}, k = 1, …, C.

The number of aggregates is ∑_{k=1}^{C}(ν_{k} + 1)(ν_{k+1} + 1)(ν_{k+2} + 1), and disaggregation is obtained by

p(z) = \frac{r_{1}(z_{1}, z_{2}, z_{3}) · . . . · r_{C}(z_{C}, z_{C+1}, z_{C+2})}{r_{2}(z_{2}, z_{3}) · . . . · r_{C}(z_{C}, z_{C+1})}, z ∈ Z^{T} .

From this one obtains the type-4 marginal probabilities

r_{k}(j, j′, j′′, j′′′) =^{△} P{Z_{k} = j, Z_{k+1} = j′, Z_{k+2} = j′′, Z_{k+3} = j′′′}
= \frac{r_{k}(j, j′, j′′)r_{k+1}(j′, j′′, j′′′)}{r_{k+1}(j′, j′′)}, k = 1, …, C − 1,

which are needed for the aggregate transition probabilities

π_{k}((j, j′, j′′), (l, l′, l′′)) =

λ_{1}Δ  if k = 1, j + 1 = l, j < ν_{k},
μ_{k−1}Δ(1 − \frac{r_{k−1}(0,j′,j′′,j′′′)}{r_{k}(j,j′,j′′)}) if k > 1, j + 1 = l, j < ν_{k},
μ_{k}Δ  if j − 1 = l, j′ + 1 = l′, j > 0, j′ < ν_{k+1},
μ_{k+1}Δ  if j′ − 1 = l′, j′′ + 1 = l′′, j′ > 0, j′′ < ν_{k+2},
μ_{k+2}Δ(1 − \frac{r_{k}(j,j′,j′′,ν_{k+3})}{r_{k}(j,j′,j′′)}) if k < C, j′′ − 1 = l′′, j′′ > 0,
μ_{k+2}Δ  if k = C, j′′ − 1 = l′′, j′′ > 0,
0  otherwise if (j, j′, j′′) ≠ (l, l′, l′′),

(j, j′, j′′), (l, l′, l′′) ∈ S_{k} × S_{k+1} × S_{k+2}, k = 1, …, C.

Clearly, the joint marginal probabilities are

r_{k}(j, j′) = \{ \begin{matrix}∑_{j′′∈S_{k+2}} r_{k}(j, j′, j′′) & if k ≤ K − 2, \\ ∑_{j′′∈S_{k}} r_{k}(j′′, j, j′) & if k > 1, \end{matrix} j ∈ S_{k}, j′ ∈ S_{k+1}, k = 1, …, K.

The blocking probabilities are obtained using (16).

2.2 Closed queueing networks

Again, we first give a Markov chain model, and then we present DA iteration with canonical and fine grain aggregation. Here the solution of nonlinear equations is needed for disaggregation.

A state z = (z_{1}, …, z_{K}) of the Markov chain gives the numbers z_{k} of jobs in node k, k = 1, …, K, and the state space is given in (15). The only transitions from state y ∈ Z^{T} into state z ∈ Z^{T} occur when a job’s service in node k is finished and this job goes to node l: z = y − 1_{k} + 1_{l} if y_{l} < ν_{l}, y_{k} > τ_{k}, and k ≠ l, with transition rate μ_{k}t_{k,l}. Thus the generator matrix Q has the elements

q_{y,z} = \{ \begin{matrix}μ_{k}t_{k,l} & if y_{k} > τ_{k}, y_{l} < ν_{l}, z = I(y − 1_{k} + 1_{l}), \\ − ∑n_{{x=1 \atop x≠y}} q_{y,x} & if y = z, \\ 0 & otherwise, \end{matrix}
y = I^{−1}(y), y, z = 1, …, n,

and P = QΔ + I with Δ = 0.99/ max_{x} |q_{x,x}| is the transition probability matrix of a Markov chain with discrete time, which has the same steady state probabilities.

Canonical aggregation

Again the number C of partitions equals the number K of nodes in the queueing network. N is the population. Aggregate Z_{k,j} is the set of all states with j jobs in node k.

The new macro probabilities R^{(t+1)} are obtained according to (8) and (9) using

φ_{k}(j, l) = \{ \begin{matrix}α_{k}^{j} & if l = j − 1 ≥ τ_{k}, \\ r_{k}(j) − β_{k}^{j} & if l = j = τ_{k}, \\ r_{k}(j) − α_{k}^{j} − β_{k}^{j} & if τ_{k} < l = j < ν_{k}, \\ r_{k}(j) − α_{k}^{j} & if l = j = ν_{k}, \\ β_{k}^{j} & if l = j + 1 ≤ ν_{k}, \\ 0 & otherwise, \end{matrix}

where

α_{k}^{j} = ∑_{{c=1 \atop c≠k}}^{K} μ_{k}Δt_{k,c}(r_{k}(j) − r_{k,c}(j, ν_{c})), β_{k}^{j} = ∑_{{c=1 \atop c≠k}}^{K} μ_{c}Δt_{c,k}(r_{k}(j) − r_{k,c}(j, τ_{c})),

j, l = τ_{k}, …, ν_{k}, k = 1, …, K.

The marginal probabilities and the normalizing constant G are efficiently calculated from the MED parameters ρ by convolution:

G = (ρ_{1} ⊗ … ⊗ ρ_{K})_{N},

r_{k}(j) = \frac{ρ_{k,j}}{G} (⨂_{{c=1 \atop c≠k}}^{K} ρ_{c})_{N−j} , j = τ_{k}, …, ν_{k}, k = 1, …, K,

r_{k,k′}(j, j′) = \frac{ρ_{k,j}ρ_{k′,j′}}{G} ( ⨂_{{c=1 \atop c≠k,k′}}^{K} ρ_{c})_{N−j−j′} , (j, j′) ∈ S_{k} × S_{k′}, \{k, k′\} ⊆ [1 : K],

where ρ_{c} =^{△} (ρ_{c,τ_{c}}, …, ρ_{c,ν_{c}}), c = 1, …, K, and ⊗ is the convolution operator (ρ_{k,j} = 0 for all j ∉ [τ_{k} : ν_{k}]). Approximate disaggregation is done according to (14).

Fine grain aggregation

Now we consider ACP aggregation for closed queueing networks with cyclic topology: The N jobs travel from node k to node k ⊕ 1 where k ⊕ 1 equals 1 if k = K and k + 1 otherwise.

Two subsequent nodes k, k ⊕ 1 may assume the states

S_{k}^{(2)} =^{△} {(z_{k}, z_{k⊕1}) | z ∈ Z^{T} }

= \{(j, j′) | max (τ_{k⊕1}, N − ∑_{{i≠k, \atop i≠k⊕1}} ν_{i} − j) ≤ j′ ≤ min (ν_{k⊕1}, N − ∑_{{i≠k, \atop i≠k⊕1}} τ_{i} − j), j ∈ S_{k}\}.

There are C = K partitions of the state space Z^{T} . The aggregate probabilities are the

joint marginal probabilities r_{k}(j, j′) =^{△} P{Z_{k} = j, Z_{k⊕1} = j′}, (j, j′) ∈ S_{k}^{(2)}, k = 1, …, K.

The state probabilities are

p(z) = ρ_{1,z_{1},z_{2}} · . . . · ρ_{K,z_{K},z_{1}}/G, z ∈ Z^{T} ,

where the MED parameters ρ_{k,j,j′}, (j, j′) ∈ S_{k}^{(2)}, k = 1, …, K, and the normalization constant G are solutions of the nonlinear equations (3) and (4).

The normalization constant G is efficiently calculated as follows, given the MED parameters. Let

σ_{k}(z_{1}, z_{k−1}, j_{k}) =^{△} \{ \begin{matrix}ρ_{K−1,z_{k−1},j_{k}}ρ_{K,j_{k},z_{1}} if k = K, \\ ∑_{z_{k}=τ_{k}}^{min(ν_{k},j_{k})} ρ_{k−1,z_{k−1},z_{k}}σ_{k+1}(z_{1}, z_{k}, j_{k} − z_{k}) if k < K, \end{matrix}

z_{1} ∈ S_{1}, j_{k} = N − z_{1} if k = 2,
z_{k−1} = τ_{k−1}, …, min(ν_{k−1}, N − z_{1}) and
j_{k} = max(0, N − z_{1} − ν_{2} − … − ν_{k−2} − z_{k−1}, τ_{k} + … + τ_{K})
, …, min(N − z_{1} − τ_{2} − … − τ_{k−2} − z_{k−1}, ν_{k} + … + ν_{K}) if k > 2,
k = 2, …, K,
σ_{1} =^{△} ∑_{z_{1}=τ_{1}}^{ν_{1}} σ_{2}(z_{1}, z_{1}, N − z_{1}).  (23)

Then G = σ_{1} holds.

The type-3 marginal probabilities,

r_{k}(j, j′, j′′) =^{△} P{Z_{k} = j, Z_{k⊕1} = j′, Z_{k⊕2} = j′′},

(j, j′, j′′) ∈ {(z_{k}, z_{k⊕1}, z_{k⊕2})|z ∈ Z^{T} }, where k⊕2 =^{△} (k⊕1)⊕1, are needed for the aggregate transition probabilities. They are calculated similarly: Replace in (23) τ_{k} and ν_{k} by j, τ_{k⊕1} and ν_{k⊕1} by j′, τ_{k⊕2} and ν_{k⊕2} by j′′ , and calculate σ_{1}; then r_{k}(j, j′, j′′) = σ_{1}/G holds.

The aggregate transition probabilities are

π_{k,k⊕1}((j, j′), (l, l′)) =

μ_{k⊖1}Δ(1 − r_{k⊖1}(τ_{k⊖1}, j, j′)/r_{k}(j, j′)) if j + 1 = l, j′ = l′, j < ν_{k},
μ_{k}Δ  if j − 1 = l, j′ + 1 = l′, j > τ_{k}, j′ < ν_{k⊕1},
μ_{k⊕1}Δ(1 − r_{k⊕1}(j, j′, ν_{k⊕2})/r_{k}(j, j′)) if j = l, j′ − 1 = l′, j′ > τ_{k⊕1},
0  otherwise if (j, j′) ≠ (l, l′),

(j, j′), (l, l′) ∈ S_{k}^{(2)}, k = 1, …, K,

where k ⊖ 1 equals K if k = 1 and k − 1 otherwise.

When DA iteration for stationary probabilities is to be done, we omit as usual the solution of the nonlinear equations (3) and (4) and apply iteration (14) instead.

The blocking probabilities are obtained using (16).

3 Numerical experiences

In order to find out what the merits of DA iteration are with respect to accuracy and computing time, we analyzed many Markov models applying this method, and compared the obtained approximations with the exact figures from the Markov chain solution. Besides we considered heuristic techniques from the literature, calculated solutions with them, and compared the accuracy of these results with ours.

Typically, an open or closed queueing network is analyzed for a variety of system parameters like number of nodes, number of jobs in the network, buffer capacity of a node, arrival rates, or service rates, using some different solution techniques like vector iteration for exact results, DA iteration with canonical aggregation and with fine grain aggregation, or methods from the literature. For each approximation technique, we report the largest observed approximation error. Furthermore, for a quite large model, we report the computing times which were consumed by vector iteration and by DA iteration with canonical and with fine grain aggregation. For a very large model which cannot be analyzed exactly, we present the computing time for DA iteration with canonical aggregation.

We present an outline of our results analyzing the following queueing network models with repetitive-service blocking:

1. Open tandem networks, arrivals at the first node. Comparison with results according to Caseau and Pujolle [3]. Their technique relies on a heuristic: the sevice rates μ_{k}, k = 1, …, K, are augmented in order to take into account repeated services which are due to blocking.

2. Open tandem networks. Comparison with results according to Labetoulle and Pujolle [10]. They use diffusion approximation for the solution. The blocking networks are open, interarrival times and service times are generally distributed. Here, the blocking policy is a bit different, a rejected job clings to the next node (repetitive-service with fixed destination). In a tandem network both blocking policies are equivalent.

3. Open networks with various topologies. Comparison with results according to Kouvatsos, Xenios and Denazis [9, 8]. Their algorithm is quite general, suited for open and closed queueing networks with repetitive-service blocking. Interarrival times and service times have general exponential distributions. The state probabilities are obtained by maximizing the entropy observing some reasonable constraints.

4. Comparison of DA iteration with coarse, fine, and very fine aggregation for open networks with tandem, split, merge, and symmetrical topology.

5. Closed, cyclic networks. Comparison with results obtained by the technique of Frein and Dallery [4]. Their decomposition method can analyze closed cyclic networks with blocking-before-service policy which is equivalent to the repetitiveservice policy in these networks.

6. Comparison of the computing times for vector iteration and DA iteration, applied to a symmetrical topology with a large number of states.

7. A tandem topology with 10^{9} system states. DA iteration consumes only a moderate amount of computing time for this large model.

The accuracy is measured using absolute errors for probabilities and utilizations and tolerances for mean job numbers and throughputs. The absolute error of any performance measure in node k is |e_{k} − a_{k}| where e_{k} is the true value of this measure for node k, and a_{k} the calculated approximation of it. The tolerance is defined as |e_{k} − a_{k}|/ ∑_{l=1}^{K} |e_{l}|, k = 1, …, K.

1. The considered queueing networks have a tandem topology with K nodes, where jobs from node k go to node k + 1, k = 1, …, K − 1. Arrivals occur only in node 1, and the jobs which leave node K leave the network.

We report the maximum absolute errors and the maximum tolerances, maximum over all nodes and all examples 1.1 and 1.2, respectively, achieved by the heuristic technique from the literature (L) and by DA iteration with canonical aggregation (DAI-canonical) and center aggregation (DAI-center):

Examples 1.1 (Tandem network) K identical nodes: λ_{1} = 1, μ_{k} = 1, ν_{k} = 3, k = 1, …, K, K = 2, …, 7. The maximum absolute errors of the probabilities of the job numbers in the nodes (macro probabilities) are 0.14 (L) and 0.007 (DAI-canonical), and for of the utilization of the servers 0.08 (L) and 0.005 (DAI-canonical), respectively. The maximum tolerances of the mean job numbers in the nodes are 0.1 (L) and 0.003 (DAI-canonical).

Examples 1.2 (Tandem network) K = 4, ν_{k} = 3, k = 1, …, K, λ_{1} = 1, μ_{1} = μ_{2} = μ_{4} = 1, μ_{3} = 2, 3, …, 10. Here we applied also DA iteration with center aggregation - center is node 3 (DAI-center). Maximum absolute errors of the macro probabilities: 0.6 (L), 0.12 (DAI-canonical), 0.02 (DAI-center), and of the utilization of the servers: 0.09 (L), 0.01 (DAI-canonical), 0.01 (DAI-center). Maximum tolerances of the mean job numbers: 0.25 (L), 0.05 (DAI-canonical), 0.01 (DAI-center).

2. Examples 2.1 The same tandem network as in examples 1.1 is investigated. The number of nodes and the capacities vary as follows: K = 1, …, 7, capacity ν = 2, 3, 4 in all nodes. The macro probabilities r_{k}(0) and r_{k}(ν), k = 1, …, K, are calculated. The maximum absolute errors, maximum over the nodes, is fairly varying among the different models (different values K) when the diffusion approximation from the literature is applied, the overall maximum is 0.57. Applying DA iteration with canonical aggregation, it is 0.013.

3. We solved some queueing netwoks with the maximum entropy algorithm from the literature (L) and with DA iteration (DAI) applying canonical or center aggregation as appropriate.

Examples 3.1 The tandem network, K = 2, …, 7 nodes, capacity ν = 3 in all nodes, and also the service rates are equal in all nodes, varying between 2 and 5.

Examples 3.2 The tandem network, K = 4, ν = 3, μ_{1} = μ_{2} = μ_{4} = 1, μ_{3} = 2, 3, …, 10.

Examples 3.3 A merge topology: K = 4 nodes, ν_{1} = 5, ν_{2} = 3, ν_{3} = 2, ν_{4} = 3, arrivals at nodes 2,3, and 4: λ_{2} = 1, λ_{3} = 0.5, λ_{4} = 0.5, transfer probabilities t_{1,0} = t_{2,1} = t_{3,1} = t_{4,1} = 1, service rates μ_{2} = μ_{3} = μ_{4} = 1, μ_{1} = 1.1, 1.2, …, 2.5.

Examples 3.4 A split topology. K = 4 nodes, ν_{1} = 5, ν_{2} = 3, ν_{3} = 2, ν_{4} = 3, arrivals at node 1, λ_{1} = 0.1, 0.2, …, 1.9, transfer probabilities t_{1,2} = 0.2, t_{1,3} = 0.3, t_{1,4} = 0.5, t_{2,0} = t_{3,0} = t_{4,0} = 1, service rates μ_{1} = 1, μ_{2} = 0.6, μ_{3} = 0.4, μ_{4} = 0.3.

Again we give the maximum absolute errors of the macro probabilities and utilizations and the maximum tolerances of the mean job numbers, maximum over the nodes and over all examples 3: Probabilities of the job numbers in the nodes: 0.23 (L), 0.02 (DAI), utilization of the servers: 0.28 (L), 0.02 (DAI), mean job numbers in the nodes: 0.21 (L), 0.02 (DAI).

All examples were solved using DA iteration with canonical aggregation except examples 3.2. Here DA iteration with canonical aggregation gives the errors 0.12, 0.02, 0.046, respectively, so center aggregation was applied in order to gain the same order of accuracy as in the other examples.

4. These numerical examples are analyzed using DA iteration with canonical and center, paired, ACP, or ACT aggregation in order to demonstrate the merits of fine grain aggregation. We report the maximum absolute errors of the macro probabilities.

Examples 4.1 The tandem network: K = 4 nodes, capacity ν = 3 in all nodes, arrivals at node 1, λ_{1} = 1, service rates μ_{1} = μ_{2} = μ_{4} = 1, μ_{3} = 1, 2, …, 20, transfer probabilities t_{1,2} = t_{2,3} = t_{3,4} = t_{4,0} = 1.

The maximum absolute error occurs when μ_{3} = 20. It is 0.13 using canonical aggregation, 0.02 using center aggregation, 0.022 using ACP aggregation, and if ACT aggregation is applied, which is even finer, it is 0.0023.

Examples 4.2 In the split network, examples 3.4, the maximum absolute error is 0.022 when canonical aggregation is used. It occurs for λ_{1} = 0.5. When center aggregation is used, it is 0.002 and occurs for λ_{1} = 0.6.

Examples 4.3 In the merge network from example 3.3, but with μ_{2} = μ_{3} = 0.5, the maximum absolute error is 0.016 when canonical aggregation is used. It occurs for μ_{1} = 2.0. When center aggregation is used, it is 0.0015 and occurs for μ_{1} = 2.2.

Examples 4.4 A symmetrical network: K = 4 nodes, capacities ν_{1} = ν_{3} = 3, ν_{2} = ν_{4} = 2, arrivals at nodes 1 and 3, λ_{1} = λ_{3} = 0.01, 0.02, …, 0.2, transfer probabilities t_{1,2} = t_{3,4} = 0.8, t_{1,3} = t_{3,1} = 0.2, t_{3,0} = t_{4,0} = 1, service rates μ_{1} = μ_{3} = 0.5, μ_{2} = μ_{4} = 0.1.

The maximum absolute error is 0.11 when canonical aggregation is used. It occurs for λ_{1} = λ_{3} = 0.08. When paired aggregation is used, K_{2} = { (1, 2), (3, 4) }, it is 0.011 and occurs also for λ_{1} = λ_{3} = 0.08.

5. We analyzed the following closed cyclic networks with the decomposition method from the literature (L) and using DA iteration with canonical aggregation (DAI-canonical) and with ACP aggregation (DAI-ACP):

Examples 5.1 K = 4 nodes, capacities ν_{1} = 5, ν_{2} = 4, ν_{3} = 5, ν_{4} = 4, service rates μ_{1} = 3, μ_{2} = 4, μ_{3} = 5, μ_{4} = 2, N=5,…,9 jobs.

Examples 5.2 K = 4 nodes, capacities ν_{1} = 2, ν_{2} = 4, ν_{3} = 3, ν_{4} = 4, service rates μ_{1} = 2, μ_{2} = 1, μ_{3} = 3, μ_{4} = 4, N=3,…,7 jobs.

Examples 5.3 K = 5 nodes, capacities ν_{1} = 2, ν_{2} = 4, ν_{3} = 3, ν_{4} = 3, ν_{5} = 4, service rates μ_{1} = 2, μ_{2} = 2, μ_{3} = 2, μ_{4} = 2, μ_{5} = 2, N=3,…,8 jobs.

Examples 5.4 K = 5 nodes, capacities ν_{1} = 2, ν_{2} = 4, ν_{3} = 3, ν_{4} = 3, ν_{5} = 4, service rates μ_{1} = 5, μ_{2} = 3, μ_{3} = 4, μ_{4} = 6, μ_{5} = 2, N=3,…,8 jobs.

The throughputs, the mean numbers of jobs in the nodes, and the macro probabilities were calculated. We give maximum tolerances and absolute errors, respectivly, maximum with respect to all nodes and all examples. Maximum tolerances of the throughputs: 0.015 (L), 0.0095 (DAI-canonical), 0.00033 (DAI-ACP), and of the mean job numbers: 0.158 (L), 0.043 (DAI-canonical), 0.0005 (DAI-ACP). Maximum absolute errors of the macro probabilities: 0.475 (L), 0.087 (DAI-canonical), 0.004 (DAI-ACP).

Example 6 This example gives an impression of the consumed computing time of DA iteration vs. vector iteration which calculates “exact” solutions. The symmetrical network of example 4.4, but with capacities ν_{1} = ν_{3} = 12, ν_{2} = ν_{4} = 8, has 13 689 states. The parameters λ_{1} = λ_{3} = 0.32, μ_{1} = μ_{3} = 2, μ_{2} = μ_{4} = 0.4 were used. On a SUN-4/330, vector iteration took 2 752 seconds, DA iteration with canonical aggregation 7 seconds, and with paired aggregation 99 seconds.

Example 7 The tandem topology, example 4.1, but with 9 nodes, μ_{3} = 1, and the capacity ν = 9 in all nodes, which has 10^{9} system states, was analyzed. On an Atari 1040ST, the computing time for DA iteration with canonical aggregation was 238 seconds. This demonstrates that Markov models with a really large state space can be treated with the DA iteration technique.

Conclusion

The DA iteration technique is well suited for the analysis of queueing networks with repetitive-service blocking, but, in fact, it is much more general, applicable to many Markov models. Its time and space demands are higher than for many special purpose heuristics from the literature, but are by orders of magnitude lower than for exact solutions of Markov chains. Hence Markov models can be treated with DA iteration whose exact analysis is impossible due to time and space complexity. Further work has to be done about the approximation error estimates.

Acknowledgements

The authors would like to express their gratitude to the unknown referees for their valuable suggestions.

References

[1] Ian F. Akyildiz and Horst von Brand. Exact solutions for open, closed, and mixed queueing networks with rejection blocking. Theoretical Computer Science, 64:203–219, 1989.

[2] B. Bärk. Aggregationsmethoden für Markovketten. 1997.

[3] P. Caseau and G. Pujolle. Throughput capacity of a sequence of queues with blocking due to finite waiting room. IEEE Transactions on Software Engineering, 5(6):631–642, 1979.

[4] Y. Frein and Y. Dallery. Analysis of Cyclic Queueing Networks with Finite Buffers and Blocking Before Service. Performance Evaluation, 10:197–210, 1989.

[5] W. K. Grassmann. Finding transient solutions in Markovian event systems through randomization. In Numerical Solution of Markov Chains, pages 357–372. Marcel Dekker, New York, Basel, Hong Kong, 1991.

[6] V. Jonas. Iterative Analyse von Blockiernetzen mit alternierender Aggregation und Disaggregation. Master’s thesis, Universität Bonn, 1992.

[7] A. M. Kagan, J. V. Linnik, and C. R. Rao. Characterization Problems in Mathematical Statistics. Wiley, New York, 1973.

[8] D. Kouvatsos and S. Denazis. Comments on and tuning to: Mem for arbitrary queueing networks with repetitive-service blocking under random routing. Research Report CS-18-91, Bradfod, 1991.

[9] D. Kouvatsos and N. P. Xenios. Mem for arbitrary queueing networks with multiple general servers and repetitive-service blocking. Performance Evaluation, 10:169–195, 1989.

[10] J. Labetoulle and G. Pujolle. Isolation method in a network of queues. IEEE Transactions on Software Engineering, 6(6):373–380, 1980.

[11] R. O. Onvural. Survey of closed queueing networks with blocking. ACM Computing Surveys, 22(2):83–121, 1990.

[12] H. G. Perros. Approximation algorithms for open queueing networks with blocking. In Takagi, editor, Stochastic Analysis of Computer and Communication, pages 451–498. North-Holland, 1990.

[13] J. Ch. Strelen. Iterative Analyse von Markov-Modellen mit alternierender Aggregation und Disaggregation. In A. Lehmann and F. Lehmann, editors, Messung, Modellierung und Bewertung von Rechensystemen IFB 286, pages 320–336. Springer, Berlin, 1991.

[14] J. Ch. Strelen. Approximate product form solutions for Markov chains. Performance Evaluation, 30:87–110, 1997.