Approximate Product Form Solutions for Markov Chains

Johann Christoph Strelen Rheinische Friedrich–Wilhelms–Universität Bonn

Abstract

An approximation technique for finite state Markov chains is proposed. The solutions are in product form and can be calculated iteratively. This iteration is equivalent to disaggregating macro probabilities, multiplying with transition probabilities, and aggregating the obtained state probabilities – we call it DA iteration. Large Markov chains with many states can be solved, so large that the state probabilities could neither be stored nor calculated in any computer.

Accurate approximations are obtained, and moreover, the accuracy can be influenced: The more aggregates are considered, the better the approximations are, if the aggregates are chosen suitably.

Practical experience is encouraging: The obtained accuracy compares favourably with other approximation methods and with exact solutions. In this paper, we discuss steady state probabilities.

Introduction

Markov chains are widely used in stochastic modeling. We are interested in models of computer, communication, and manufacturing systems. There are discrete time and continuous time Markov chains. The latter can be transformed into discrete time Markov chains which behave equivalently. This is why our technique, which is designed for discrete time Markov chains, is useful for both types.

The goal of considering Markov chain models is to obtain performance measures like throughput of jobs, utilization and availability of system components, and response times. Usually, throughput, utilization, availability, and in some cases response times can be derived directly from state probabilities, but sometimes this is difficult for response times ([2],[5]). Steady state probabilities are solutions of typically large and sparse linear systems of equations.

The most important, well-known, and widely applied modeling paradigms in this field, namely queueing systems, queueing networks, and generalized stochastic Petri nets (GSPN), can be considered as abstract representations of Markov chains, if the memoryless condition is satisfied – this can be assumed if all probability distributions are of Coxian type.

0This work was supported in part by the Deutsche Forschungsgemeinschaft

Unfortunately, the number of Markov chain states grows quickly as the complexity of the model increases. Therefore much effort has been made on the development of algorithms for large linear systems of equations, like power, Gauss-Seidel, and SOR algorithms; an overview is given in [11]. Nowadays, systems with some hundreds of thousands of states can be solved. A special class of algorithms is based on the principle of iterative aggregation and disaggregation, a well-established solution technique for Markov chains, see for example the survey of Schweitzer [12].

Another problem are different time scales which stem from rare events in the modeled systems. This may lead to very long iteration times when iterative algorithms are applied. Simultaneous iterative techniques ([6], [14], [15], [16]) cope with this.

Our 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; we call them disaggregation aggregation product form (DAΠ) solutions. 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.

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.

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.

In this paper, we consider steady state solutions. In section 1, DAΠ solutions and DA iteration are discussed. In section 2, we show that there exists always a product form approximation for any discrete probability distribution, which is as accurate as desired, with only a small number of factors. Practical aspects are treated in section 3. As an illustrative example, the techniques are applied to queueing networks with blocking.

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.

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 is an aggregate.

If for all state tuples z = (z_{1}, …, z_{C}) ∈ [0 : N_{1}] × … × [0 : N_{C}] the intersections Z_{1,z_{1}} ∩ … ∩ Z_{C,z_{C}} contain at most one element, the aggregation is said to completely identify the states. If the intersection is empty, the state tuple z is termed infeasible, otherwise 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.

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) for z = I(z), z ∈ Z^{T}, holds.

The n-vector p^{(t)} = [p^{(t)}(1), …, p^{(t)}(n)] contains the state probabilities, 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 (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.

Sometimes joint marginal probabilities

r_{k_{1},k_{2}}(j_{1}, j_{2}) =^{△} P{Z ∈ Z_{k_{1},j_{1}}, Z ∈ Z_{k_{2},j_{2}}} = ∑_{z∈Z_{k_{1},j_{1}}∩Z_{k_{2},j_{2}}} p(z),

j_{i} = 0, …, N_{k_{i}}, {k_{1}, k_{2}} ⊆ [1 : C],

are used. Again we collect them into tuples r_{k_{1},k_{2}}.

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 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.

A special aggregation is defined by the structure of the model as follows. Consider a system which consists of K subsystems, e.g. a queueing network with K nodes. If each aggregate is defined by the state z_{k} of one 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.

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 [8, 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 parameters ρ_{k,j}, ρ for short.

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

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

In general, the nonlinear equations (3) must be solved in order to calculate all MED parameters ρ_{k,j} and G of the MED π, given the marginal probabilities R. This is easily accomplished in some special cases. For example, if the 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 [16]). Here ρ_{k,j} = r_{k}(j), j = 0, …, N_{k}, k = 1, …, C, and G = 1 hold. This is generalized in section 3 for the case where Z^{T} consists of the union of cartesian products. Another simple case 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 can 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)}.

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 (we tried quite a lot). In section 3, these matters are treated.

More precisely, the conditions are fulfilled if the distribution functions φ_{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),

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

holds. (8) or (11) are used to calculate the new macro probabilities R^{(t+1)} from the old ones, R^{(t)}. In this way the function T is evaluated.

The nonlinear equations for the MED parameters ρ, given the macro probabilities R, can be solved iteratively as follows ([7]), if no closed form solution is available:

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

Here ρ^{[0]} 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)} = \frac{r_{k}^{(t+1)}(j)}{r_{k}^{(t)}(j)}ρ_{k,j}^{(t)}  (14)

within each DA iteration step (6).

We summarize, stating the

DA Iteration Algorithm

Input: Reasonable MED parameters ρ^{(1)} – e.g. 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 using (5) or (28), 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)} and the MED parameters ρ^{(t+1)} and G^{(t+1)}. □

Now we are going to prove the existence of DAΠ solutions. This is accomplished applying the following theorem, see e.g. [10, 6.3.2].

Brouwer Fixed-Point Theorem Let T′ : P̄_{n} ⊂ IR^{n} → IR^{n} be a continuous operator on the compact, convex set P̄_{n}, and suppose that T′(P̄_{n}) ⊂ P̄_{n}. Then T′ has a fixed point in P̄_{n}.

DA Product Form Existence Theorem The equation

R = T (R)

has a fixed point.

Proof Instead of T , we consider the operator

T′ : P̄_{n} → P̄_{n}, p ↦ D(A(p)) · P.  (15)

P̄_{n} = {p′ ∈ [0, 1]^{n} | p′_{1} + … + p′_{n} = 1} ⊂ IR^{n} is the set of all probability distributions over [1 : n]. Obviously, P̄_{n} is bounded and closed, therefore compact, and convex.

Now we prove that T′ is defined and continuous. In the sequel, we only consider marginal probabilities R, for which a probability distribution p̃ exists such that R = A(p̃) holds. Therefore we may assume that always a MED exists with the constraints that the marginal probabilities R are assumed.

Case 1 Z^{T} is rectangular. Firstly, let p ∈ P̄_{n}, R =^{△} A(p). Clearly, R depends continuously on p. Secondly, let

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

p′ is a probability distribution, p′ ∈ P̄_{n}, which depends continuously on R. In the next paragraph, we show, that (16) is equivalent to p′ = D(R). Thirdly, p′P is a continuous mapping. Consequently, T′ (15) is defined in this case, and it is continuous.

If all marginal probabilities r_{k}(z_{k}) are positive, p′ is a MED, according to the product form theorem. If there are marginal probabilities r_{k}(z_{k}) which are zero, p′ is also a MED, as can be seen as follows: Let

Z̃ = Z^{T} \ {z ∈ Z^{T} | ∃k, j : z ∈ Z_{k,j} where r_{k}(j) = 0}.

According to the product form theorem,

p̃(z) = p′(z), z ∈ Z̃,  (17)

is the MED, given the positive marginal probabilities as constraints. For each state z ∈ Z_{k,j} where r_{k}(j) = 0, p′(z) = 0 must hold if the constraints R = A(p′) are given. This is fulfilled by p′ in (16). Therefore p′ maximizes the entropy and is in accordance with the constraints.

Case 2 Z^{T} is not rectangular, but a proper subset of [0 : N_{1}] × … × [0 : N_{C}]. Therefore Z^{R} =^{△} [0 : N_{1}] × … × [0 : N_{C}] \ Z^{T} is not empty. Here, we use the subset independence feature of maximum entropy distributions, which is stated in the Shore and Johnson paper [13]. For our purposes, this can be formulated as follows. Let p̃ denote a probability distribution over Z^{R} where all p̃(z) > 0, z ∈ Z^{R}, R′′′ =^{△} A(p̃), w a probability, 0 < w < 1, and

R′′ = wR + (1 − w)R′′′.

Clearly, R′′ = A(wp + (1 − w)p̃) holds. Therefore there is a MED p′′ of the form (16) obeying the constraints R′′ = A(p′′). According to the subset independence feature,

p′′ = wp′ + (1 − w)p′′′

holds, where p′ is the unique MED over Z^{T} with the constraints R = A(p′), and p′′′ the unique MED over Z^{R} with the constraints R′′′ = A(p′′′).

Keeping R′′′ and w fixed, R′′ depends continuously on R. According to (16), applied to R′′ and p′′, p′′ depends continuously on R′′. Observing

p′ = (p′′ − (1 − w)p′′′)/w,

p′ is defined, and depends continuously on p′′. Together, p′ depends continuously on R.

In each of the two cases, we saw that T′ is defined. Moreover, R = A(p), p′ = D(R), and p′P are continuous in each case. Therefore, T′ is continuous. Furtermore, for each p ∈ P̄_{n}, T′(p) is a probability distribution, i.e. T′(p) ∈ P̄_{n}.

So far we have proved that T′ and P̄_{n} fulfil all assumptions of the Brouwer fixed-point theorem. As a consequence, T′ has a fixed point in P̄_{n}, say p.

Let R =^{△} A(p). In turn, T (R) = A(D(R) · P) = A(D(A(p)) · P) = A(T′(p)) =

A(p) = R holds. Clearly, R is a fixed point of T . This concludes the proof. □

Unfortunately, we have not yet found error estimates or a proof of convergence of the iteration. But 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 could be substantially ameliorated if more and suitably chosen aggregates are considered. These matters are discussed in section 3.

2 Motivation

In this section, we justify why it is worth applying DA iteration to Markov models.

First, we mention favourable experience with DAΠ approximations for steady state solutions of Markov models: We calculated DAΠ approximations for many models using DA iteration. Usually, we tried different aggregations for the same model, coarse grain with a small number of aggregates, which use a small amount of computing time, and fine grain with more aggregates, which need more computing time and storage space, but produce on the other hand even more accurate approximations. In section 3, these topics are pointed out.

In [17] and [16], DA iterative solutions for polling systems are treated, in [16] for queueing networks with a fork-join subnet, in [9] and [4] for some generalized stochastic Petri nets, and in [18] for queueing networks with repetitive service blocking.

The obtained accuracy compares favourably with results calculated using other methods from the literature, and with exact results.

We describe the blocking networks in some detail:

Example 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 the service of a job in node k is finished, the job tries to go to node l with 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.

In section 3, we present the analysis of these blocking networks using DA iteration.

□

Secondly, we look on general product form approximations π(z) = ρ_{1,z_{1}} · . . . · ρ_{C,z_{C}}ρ_{0}, z ∈ Z^{T}, for probability distributions q. We will show that there is always such an approximation with a small number of factors (the MED parameters ρ_{k,z_{k}}), given any degree of accuracy, namely an ϵ′ such that ||π −q||_{1} ≤ ϵ′ holds. The proof is constructive: The aggregation of the product form approximation π and the MED parameters are stated. This is accomplished directly, not using marginal probabilities as restrictions when maximizing the entropy.

For the treatment of general product form approximations π we do not require that the aggregation identifies the states completely. Thus the mapping I^{−1} : Z → Z^{T} is not always one-to-one: Given a state z ∈ Z, and letting z = (z_{1}, …, z_{C}) = I^{−1}(z) denote the state tuple belonging to z, z ∈ Z_{k,z_{k}} holds for all k = 1, …, C, but there may exist another state z′ ∈ Z with the same tuple, z = I^{−1}(z′).

Here we restrict ourselves to partitions consisting of two aggregates each, Z_{k} = {Z_{k,0}, Z_{k,1}}, k = 1, …, C.

The following theorem states the existence of a general product form approximation π, which approximates a given probability distribution q over Z accurately. In the first part of the proof, the aggregation is provided, and the MED parameters ρ_{k,j}, k = 1, …, C, j = 0, 1, are given.

Accurate Product Form Theorem Given a probability distribution q = (q(1), …, q(n)) and a degree ϵ of accuracy, 0 < ϵ ≤ 0.01, there exists an

approximating product form distribution π = (π(1), …, π(n)) with

C ≤ Ĉ =^{△} ⌈ln (\frac{ln(n/ϵ)}{− ln(1 − 2ϵ)} + 1) / ln 2⌉  (18)

partitions, each of them with two aggregates. The accuracy is

||q − π||_{1} < 4.2ϵ,  (19)

or, more precisely,

|q(z) − π(z)| < \frac{1.1}{n}ϵ if q(z) ≤ μ,
\frac{|q(z) − π(z)|}{q(z)} < 3.1ϵ if μ ≤ q(z) ≤ 1,
|q(z) − π(z)| < 2.1ϵ if 1 − ϵ ≤ q(z) ≤ 1,
z = 1, …, n,  (20)

where

μ =^{△} (1 − 2ϵ)^{2_{C}−1}.  (21)

Proof The proof is constructive: We give the product form distribution and the aggregation belonging to it, and afterwards we will prove the quoted statements about it.

The aggregation consists of C partitions Z_{k} = {Z_{k,0}, Z_{k,1}}, each with two aggregates. Let

C = Ĉ (this is preliminary, perhaps C will decrease later on),
b = 1 − 2ϵ,
ρ_{k,0} = 1, ρ_{k,1} = b^{2^{k−1}}, k = 1, …, Ĉ,
l_{z} = ⌊ln q(z)/ ln b + 0.5⌋,  (22)

π̃_{z} = \{ \begin{matrix}μ, if & q(z) ≤ μ, \\ b^{l_{z}}, if & μ < q(z) ≤ 1 − ϵ, \\ 1, if & q(z) < 1 − ϵ, \end{matrix}  (23)
ρ_{0} = 1/(π̃_{1} + … + π̃_{n}),
π(z) = ρ_{0}π̃_{z}, z ∈ Z.  (24)

The aggregation is obtained as follows. From μ < q(z) ≤ 1 − 2ϵ, 1 ≤ l_{z} ≤ 2_{C} − 1 can easily be shown. Therefore we may represent the integers l_{z} as binary numbers,

l_{z} = l_{z,1}2^{0} + l_{z,2}2^{1} + … + l_{z,_{C}}2^{_{C}−1},

which defines the figures l_{z,k} ∈ {0, 1}, k = 1, …, Ĉ, uniquely. Thus we have by (23)

π̃_{z} = μ = ρ_{1,1} · . . . · ρ_{_{C},1}  if q(z) ≤ μ,
π̃_{z} = b^{l_{z,1}·2^{0}}b^{l_{z,2}·2^{1}} · . . . · b^{l_{z,_{C}}·2^{_{C}−1}}
= ρ_{1,l_{z,1}} · . . . · ρ_{_{C},l_{z,_{C}}}  if μ < q(z) ≤ 1 − ϵ,
π̃_{z} = 1 = ρ_{1,0} · . . . · ρ_{_{C},0}  if 1 − ϵ < q(z).  (25)

Clearly,

π̃_{z} = ρ_{1,l_{z,1}} · . . . · ρ_{_{C},l_{z,_{C}}}

where l_{z,k} =^{△} 1 if q(z) < μ, and l_{z,k} =^{△} 0 if 1 − ϵ < q(z), k = 1, …, Ĉ, is in accordance with (25).

From this representation, a preliminary aggregation is obtained: A state z ∈ Z belongs to the aggregate Z_{k,j} iff in (25) ρ_{k,j} is a factor of π̃_{z}, i.e.

Z_{k,j} = {z ∈ Z | l_{z,k} = j}, j = 0, 1, k = 1, …, Ĉ.

If an aggregate Z_{k,j} is empty, then ρ_{0} is replaced by ρ_{0}ρ_{k,1−j}, the partition Z_{k} is discarded, the numbering of the remaining partitions is adjusted to 1, …, C − 1, and C is reduced by 1.

If there are two partitions Z_{k} and Z_{k′}, for which there are numbers j, j′ ∈ {0, 1} such that Z_{k,j} = Z_{k′,j′} holds, then we replace ρ_{k,j} by ρ_{k′,j′}ρ_{k,j} and ρ_{k,1−j} by ρ_{k′,1−j′}ρ_{k,1−j}, and discard the partition Z_{k′} like just explained.

This was the more interesting part of the proof, namely the definition of the approximating product form distribution and of the aggregation belonging to it.

The second part is to prove the accuracy (20), from which (19) follows immediately. The arguments are straightforward, lengthy, and consist in applying some Taylor series with remainders and simple estimates. This second part of the proof is given in the appendix. □

Remark The aggregation has C ≤ Ĉ partitions, each of them with two aggregates. This C is also the number of factors in the approximating product form. There are C degrees of freedom, namely the MED parameters ρ_{k,1}, k = 1, …, C (the parameters ρ_{0} are unity, and the normalization constant is determined by the MED parameters). The number Ĉ is O(log log n) for ϵ fixed, and O(log 1/ϵ) for n fixed. For example, if n varies between 1000 and 10^{27}, and ϵ between 0.01 and 10^{−5}, Ĉ is not greater than 22.

As a consequence of the theorem, for any probability distribution q with n−1 degrees of freedom – n might be large – there exists an accurate product form approximation π with a small number C of degrees of freedom. This is the important feature for our motivation: It is worth searching for product form approximations.

3 Practical aspects

If a Markov model is to be analyzed with DA iteration, the aggregation, i.e. the partitions of the state space, must be defined. Moreover, in the DA iteration algorithm, given the MED parameters ρ_{k,j}, the normalization constant G, the marginal probabilities r_{k}(j), and the aggregate transition probabilities π_{k}(j, l) must be calculated in an efficient manner. The defining equations of these quantities, (4), (3), and (7), (10) are composed of too many summands.

For the definition of the aggregation, we consider the modeled system as being composed of components k, k = 1, …, K. Each system component k may assume a state z_{k} of its own state space S_{k}, e.g. S_{k} = [0 : N_{k}], N_{k} ∈ IN. Consequently, the state of the whole system can be described by K-tuples z = (z_{1}, …, z_{K}) ∈ S_{1} × … ×S_{K}. Clearly, the state space is a subset of this cartesian product, Z^{T} ⊂ S_{1} × … × S_{K}.

Example In the queueing network with blocking, the system’s components are the nodes. Their states are the numbers z_{k} of jobs in them, k = 1, …, K. Therefore, a system’s state is given by the K−tuple z = (z_{1}, …, z_{K}). In an open queueing network,

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

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

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

where

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

A state τ_{k} in node k occurs if in no other node space is left. □ In this section, we consider only aggregations which identify the states completely. In section 1, canonical aggregation is defined as follows: There are C = K partitions

Z_{k}, k = 1, …, K, and the aggregates of each partition Z_{k} are defined by the states of the system component k:

Z_{k,j} = {z ∈ Z | z_{k} = j}, j ∈ S_{k}, k = 1, …, K

(If z and z occur in the same formula, both denotations correspond: z = I^{−1}(z)).

Example In the example blocking networks, each aggregate Z_{k,j} consists of all system states where j jobs are in node k. □

If more aggregates are applied, we call the aggregation fine. This can be accomplished as follows: A partition is defined by two or more system components, e.g. the partition Z_{(k,k′)}:

Z_{(k,k′),(j,j′)} = {z ∈ Z | z_{k} = j, z_{k′} = j′}, j ∈ S_{k}, j′ ∈ S_{k′}, k ≠ k′, k, k′ ∈ [1 : K]. (27)

Example In our blocking networks, we may define aggregates according to the state of two nodes. Thus the dependencies between two such nodes are completely taken into account. In general, not all partitions Z_{(k_{1},k_{2})}, (k_{1}, k_{2}) ∈ [1 : K]^{2}, are considered, but only for a subset K_{2} ⊂ [1 : K]^{2} of node pairs.

Sometimes the dependencies between a certain node k^{∗} ∈ [1 : K] and the other nodes are of particular importance, e.g. 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^{∗}},

and the partition Z_{k^{∗}}. Node k^{∗} is called the center. In an open network, the number of aggregates is

(ν_{k^{∗}} + 1) (1 + ∑_{{c=1 \atop c≠k^{∗}}}^{K} (ν_{c} + 1)) .

Sometimes the states of distinct pairs (k_{1}, k_{2}) of nodes depend strongly. Under such circumstances one can put these pairs into K_{2} and consider the according partitions Z_{(k_{1},k_{2})}, (k_{1}, k_{2}) ∈ K_{2}. If K is even, and if 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. In an open network, the number of aggregates is

∑_{(k_{1},k_{2})∈K_{2}}(ν_{k_{1}} + 1)(ν_{k_{2}} + 1).

□

In order to obtain a fine aggregation, also a virtual system component can be used. Its states s ∈ S are also a function of the system’s states, σ : Z → S, like the states of real components. Useful for the definition of fine grain aggregations are, for example, the number of jobs in an open queueing network, i.e. s = z_{1} + … + z_{K}, z ∈ Z^{T}, or the set of empty/non-empty system components, i.e.

S = {s ∈ {0, 1}^{K} | s_{k} = sign(z_{k}), k = 1, …, K, z ∈ Z^{T} }.

In general,

S = {σ(z) | z ∈ Z}

holds.

Applying virtual system components, the aggregates may be defined by

Z_{k,(j,s)} = {z ∈ Z | z_{k} = j, σ(z) = s, }, j ∈ S_{k}, s ∈ S, k = 1, …, K.

Another partition Z_{C}, C = K + 1, is defined solely the virtual component. Its aggregates are

Z_{C,s} = {z ∈ Z | σ(z) = s}, s ∈ S.

In this case, or similarly if the system’s component C represented by the component state space S is a real one, the disaggregation is sometimes quite simple. We mean the case where the state space is the union of some disjoint cartesian products as follows:

Z = ⋃_{s∈S} Z_{C,s} with Z_{C,s} = S′_{1}(s) × … × S′_{C−1}(s) × {s}

where the sets S′_{k}(s) ⊂ S_{k} may depend on s ∈ S.

Assume the marginal probabilities

r_{k}(j, s) = P{Z ∈ Z_{k,(j,s)}}, j ∈ S_{k}, s ∈ S, k = 1, …, C − 1,

and

r_{C}(s) = P{Z ∈ Z_{C,s}}, s ∈ S,

are given. Consider a given value s ∈ S. Because of the rectangular form of the subspace Z_{C,s},

P{Z = z | Z ∈ Z_{C,s}}

= P{Z_{1} = z_{1} | Z ∈ Z_{C,s}} · . . . · P{Z_{C−1} = z_{C−1} | Z ∈ Z_{C,s}}
= \frac{r_{1}(z_{1}, s)}{r_{C}(s)} · . . . · \frac{r_{C−1}(z_{C−1}, s)}{r_{C}(s)}, z ∈ Z_{C,s}

holds. From this, for all z ∈ Z_{C,s} and all s ∈ S, in turn

P{Z = z} = P{Z = z, Z ∈ Z_{C,s}}
= P{Z = z | Z ∈ Z_{C,s}} · P{Z ∈ Z_{C,s}}
= \frac{r_{1}(z_{1}, s)}{r_{C}(s)} · . . . · \frac{r_{C−1}(z_{C−1}, s)}{r_{C}(s)} · r_{C}(s), s ∈ S,

hold. Thus

ρ_{k,(j,s)} = \frac{r_{k}(j, s)}{r_{C}(s)}, j ∈ S′_{k}(s), k = 1, …, C − 1,
ρ_{C,s} = r_{C}(s), s ∈ S,  (28)

are the MED parameters.

Example In our blocking network, if center aggregation is applied, let k^{∗} = K = C, i.e. s ∈ S are the states of the center. □

In general, the aggregation cannot be chosen as we did in the proof of the accurate product form theorem, as

– the exact solution is not known and

– the obtained aggregation could be of such a kind that the normalization constant G or the aggregate transition probabilities π_{k}(j, l) cannot be calculated efficiently.

Instead, nearly dependent and independent components of the system are identified. Dependent components should be combined, in order to define together a partition as indicated in (27). Vice versa, independent components may define different partitions: If the state space is rectangular, and the components are completely independent, then the according DAΠ solution is exact and the transition probability matrix of the Markov chain is a Kronecker product, see e.g. [3]. Bärk [1] does some research on these topics.

If the state space is rectangular, Z^{T} = S_{1}× … ×S_{C}, e.g. in an open queueing network, the normalization constant

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

is efficiently calculated by

G = ( ∑_{j∈S_{1}} ρ_{1,j}) · . . . · ( ∑_{j∈S_{C}} ρ_{C,j}) .  (29)

This technique can easily be generalized if the state space is the union of some disjoint cartesian products.

We call a state space

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

where S_{k} = [u_{k} : o_{k}], a closed form state space, because closed queueing networks have a state space of this kind. In this case, the normalization constant is efficiently calculated using convolution. Let ρ_{k}, k = 1, …, C, denote the vectors (ρ_{u_{k}}, …, ρ_{o_{k}}) of DME parameters. The normalization constant G is given by the Nth component of the convolution vector of all vectors ρ_{k}, k = 1, …, C:

G = (ρ_{1} ⊗ … ⊗ ρ_{C})_{N}.  (31)

There are more general state spaces, for which both simplification techniques are applied.

The efficient calculation of a marginal probability r_{k}(j) (3) can be done accordingly: In formula (4) for the normalization constant G, ρ_{k,l} is replaced for all l ≠ j by zero, and the sum is divided by G.

More general, probabilities P{Z ∈ X } = ∑_{z∈X} ρ_{1,z_{1}} · . . . · ρ_{C,z_{C}}/G where Z = (Z_{1}, …, Z_{C}), X = X_{1} × … × X_{C}, X_{k} ⊂ S_{k}, k = 1, …, C, can efficiently be calculated as follows:

If the state space is rectangular,

P\{Z ∈ X \} = \frac{1}{G} ∏_{k=1}^{C} P\{Z_{k} ∈ X_{k}\}  (32)

holds. If the state space is closed form (30),

P{Z ∈ X ∩ Z^{T} } = (x_{1} ⊗ … ⊗ x_{C})_{N}/G  (33)

holds where x_{k} = (x_{k,0}, …, x_{k,N_{k}}),

x_{k,j} =^{△} \{ \begin{matrix}ρ_{k,j}, & if j ∈ X_{k}, \\ 0 & otherwise, \end{matrix} j ∈ S_{k}, k = 1, …, C.

The proof of the convolution formula (33) is straightforward.

Now we turn to the problem of efficiently calculating the aggregate transition probabilities. In general, this task needs O(n^{2}) operations. This is why we restrict the considered class of state transition matrices.

We assume that the Markov chain transitions from state y into state z, y ≠ z, y, z ∈ Z^{T}, can be classified as follows: Each state transition belongs to one and only one of the transition classes, which are numbered 1, …, ν̄.

A transition y into z of transition class ν may occur only if y ∈ U_{ν}, where U_{ν} denotes the source states set of this class. The destination state function u_{ν} : Z^{T} → Z^{T} defines the new states, i.e. a transition of class ν leads from state y ∈ U_{ν} into state u_{ν}(y). u_{ν} consists of component functions u′_{ν,k},

u_{ν} = (u′_{ν,1}(y), …, u′_{ν,C}(y)).  (34)

Transitions of class ν occur with probability α_{ν}.

In order to reduce the complexity of the calculation of the aggregate transition probabilities, we restrict the sets of source states to

U_{ν} = U_{ν,1} × … × U_{ν,C} ⋂ Z^{T}  (35)

where

U_{ν,k} ⊂ S_{k}, k = 1, …, C, ν = 1, …, ν̄,

and we assume, that the destination state functions u′_{ν,k} depend only on y_{k}. In this case, there are functions u_{ν,k} : U_{ν,k} → S_{k}, for which u_{ν,k}(y_{k}) = u′_{ν,k}(y), y ∈ U_{ν}, holds. Therefore we may consider

u_{ν}(y) = (u_{ν,1}(y_{1}), …, u_{ν,C}(y_{C}))  (36)

instead of (34).

The introduction of transition classes not only helps to obtain the lower complexity, but also serves to structure the design of the Markov models.

Example In an open queueing network, the system’s state may change due to an arrival at node k from outside. The source state set is

U_{ν} = {z ∈ Z^{T} | z_{k} < ν_{k}}.

That means, there is space left in node k for an arriving job. The destination state function is given by

u_{ν,k}(z_{k}) = z_{k} + 1,
u_{ν,i}(z_{i}) = z_{i},  if i ≠ k.

The transition probability is

α_{ν} = λ_{k}Δ.

Δ is the randomization constant: Initially, the Markov chain is continuous in time. Let Q denote its generator matrix. We consider the discrete time Markov chain with the transition probability matrix P = QΔ + I instead, where Δ = 0.99/ max_{x} |q_{x,x}|. It has the same steady state probabilities as the continuous time Markov chain. By the factor 0.99, numerical problems are avoided, which may occur due to rounding errors.

Other state transitions take place when a job of node k leaves the network. The according transition class ν is defined as follows. The source state set is

U_{ν} = {z ∈ Z^{T} | z_{k} > 0},

the destination state function

u_{ν,k}(z_{k}) = z_{k} − 1,
u_{ν,i}(z_{i}) = z_{i},  if i ≠ k,

and the transition probability

α_{ν} = μ_{k}(1 − t_{k,1} − … − t_{k,K})Δ.

The third kind of state transitions occur when a job leaves node k for node i, i ≠ k. In node k, there must be at least one job, and in node i must be buffer space available. Therefore, the source state set is

U_{ν} = {z ∈ Z^{T} | z_{k} > 0, z_{i} < ν_{i}}.

The states of the nodes k and i change,

u_{ν,k}(z_{k}) = z_{k} − 1,
u_{ν,i}(z_{i}) = z_{i} + 1,
u_{ν,m}(z_{m}) = z_{m},  if m ≠ i, k,

and the transition probability is

α_{ν} = μ_{k}t_{k,i}Δ.
□

Applying the transition classes, the state transition probabilities of the Markov chain can be expressed as follows:

p(y, z) = ∑_{ν=1}^{ν̄} I(y ∈ U_{ν})I(u_{ν}(y) = z)α_{ν}, y ≠ z, y, z ∈ Z^{T} ,
p(y, y) = 1 − ∑_{ν=1}^{ν̄} I(y ∈ U_{ν})α_{ν}, y ∈ Z^{T} ,  (37)

where I : {true, false} → {0, 1} denotes the indicator function.

From the definition (7), (10) of the aggregation transition probabilities, one obtains for j ≠ l, using the expressions (37) for the transition probabilities, the sets U_{ν} of source states (35), and the destination state functions u_{ν} (36), in turn

π_{k}(j, l)r_{k}(j)
= ∑ ; y∈Z^{T}_{y_{k}=j} ∑ ; z∈Z^{T}_{z_{k}=l} p(y, z)p(y)

= ∑ ; y∈Z^{T}_{y_{k}=j} ∑ ; z∈Z^{T}_{z_{k}=l} ∑_{ν=1}^{ν̄} I(y ∈ U_{ν})I(u_{ν}(y) = z)α_{ν}p(y)

= ∑_{ν=1}^{ν̄} α_{ν} ∑ ; y∈Uν_{y_{k}=j} I(u_{ν,k}(y_{k}) = l)p(y)

= ∑_{ν=1}^{ν̄} \{ \begin{matrix}α_{ν}P\{Z ∈ U_{ν,1} × … × U_{ν,k−1} × \{j\} × U_{ν,k+1} × … × U_{ν,C}\}, \\ if j ∈ U_{ν,k}, u_{ν,k}(j) = l, \\ 0,  otherwise. \end{matrix}  (38)

The probability in formula (38) can be evaluated efficiently using the formulae (32) or (33), respectively.

The aggregate transition probabilities π_{k}(j, j) are efficiently calculated using (12). Example Here, we consider only closed queueing networks with blocking. The

only transitions, which occur, are of the third kind: a job changes from a node to another. Observing the actual source state sets, destination state functions, and transition probabilities, one obtains the aggregat transition probabilities, using (38),

π_{k}(j, l)r_{k}(j) = \{ \begin{matrix}∑_{k≠i} μ_{k}t_{k,i}ΔP\{Z_{k} = j, Z_{i} < ν_{i}\} if & j > τ_{k}, l = j − 1, \\ ∑_{k≠i} μ_{i}t_{i,k}ΔP\{Z_{k} = j, Z_{i} > τ_{i}\} if & j < ν_{k}, l = j + 1. \end{matrix}

Clearly, P{Z_{k} = j, Z_{i} < ν_{i}} equals P{Z_{k} = j} − P{Z_{k} = j, Z_{i} = ν_{i}} (observe that the last probability is 0 if (j + ν_{i}) ∉ [τ_{k} + τ_{i} : N]).

If |j − l| > 1, the aggregate transition probability is 0. □ At last we mention some numerical results for open queueing networks with repetitive

service blocking.

Example A 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 of the macro probabilities is 0.13 with canonical aggregation, and 0.02 with center aggregation. It occurs when μ_{3} = 20. □

Example 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.

The maximum absolute error of the macro probabilities 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. □

Example 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} = 0.5, μ_{4} = 1, μ_{1} = 1.1, 1.2, …, 2.5.

The maximum absolute error of the macro probabilities 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. □

Example 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 also occurs for λ_{1} = λ_{3} = 0.08. □

Example The tandem topology was analyzed, but here with 9 nodes, μ_{3} = 1, and the capacity ν = 9 in all nodes, which has 10^{9} system states. 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

Product form approximations for the steady state probabilities of Markov chains and a technique for calculating them are proposed, which are suitable for huge state spaces. Unfortunately, we have not yet found error estimates or a proof of convergence of the iteration. But 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, by choosing aggregations with more aggregates – fine grain aggregation – the accuracy can be substantially ameliorated.

The just mentioned open questions, different time scales, the development of a tool, and time-dependent probabilities are topics of future research. The latter can be obtained via DA iteration: If the Markov chain is discrete in time, the DA solutions p^{(t)}, t = 1, 2, …, are approximate time-dependent state probabilities. If the Markov chain is continuous in time, randomization can be applied onto the DA solutions in order to obtain approximate time-dependent solutions.

Appendix

Proof of the accurate product form theorem, part 2.

We begin estimating the number μ = b^{2_{C}−1}: From

Ĉ ≥ ln (\frac{ln(n/ϵ)}{− ln b} + 1) / ln 2,

according to (18), we obtain

b^{2_{C}−1} ≤ ϵ/n,

and using the definition of μ,

μ ≤ ϵ/n.  (39)

Case q(z) ≤ μ: From π̃_{z} = μ, |π̃_{z} − q(z)| ≤ μ follows, and using (39), one gets

|π̃_{z} − q(z)| ≤ ϵ/n.  (40)

Case μ < q(z) ≤ 1 − ϵ: Let

r_{z} =^{△} \frac{q(z) − π̃_{z}}{q(z)} = 1 − π̃_{z}/q(z) = 1 − b^{l_{z}−c_{z}},  (41)

where c_{z} =^{△} ln q(z)/ ln b, i.e. q(z) = b^{c_{z}}. l_{z} = ⌊c_{z} + 0.5⌋ follows using (22), and in turn

c_{z} − 0.5 < l_{z} ≤ c_{z} + 0.5, −0.5 < l_{z} − c_{z} ≤ 0.5, b^{−0.5} > b^{l_{z}−c_{z}} ≥ b^{0.5},

and using (41)

1 − b^{−0.5} < r_{z} ≤ 1 − b^{0.5}.  (42)

Considering the Taylor series

b^{−0.5} = (1 − 2ϵ)^{−0.5} = 1 + ϵ + \frac{3}{2}(1 − θ)^{−5/2}ϵ^{2} where 0 ≤ θ ≤ 2ϵ,

and observing 0 ≤ ϵ ≤ 0.1,

b^{−0.5} ≤ 1 + 1.0157 … ϵ

holds, and accordingly

b^{0.5} = (1 − 2ϵ)^{0.5} = 1 − ϵ − \frac{1}{2}(1 − θ)^{−3/2}ϵ^{2} where 0 ≤ θ ≤ 2ϵ,

b^{0.5} ≥ 1 − 1.0051 … ϵ.

Applying this to (42),

−1.0157 … ϵ < r_{z} ≤ 1.0051 … ϵ,
|r_{z}| < r̂ =^{△} 1.016ϵ,  (43)

and using (41),

|q(z) − π̃_{z}| < r̂q(z)  (44)

follow.

Case 1 − ϵ < q(z): Here π̃_{z} is unity, and therefore

|π̃_{z} − q(z)| < ϵ  (45)

holds.

For calculating the approximation errors, we need bounds for the normalization constant ρ_{0} = 1/(π̃_{1} + … + π̃_{n}):

Case where there is no state z ∈ Z with q(z) > 1 − ϵ:

ρ_{0}^{−1} = ∑_{{z∈Z \atop q(z)≤μ}} π̃_{z} + ∑_{{z∈Z \atop q(z)>μ}} q(z)(1 − r_{z}) (using (41))
= ∑_{z∈Z} q(z) + ∑_{{z∈Z \atop q(z)≤μ}} (π̃_{z} − q(z)) − ∑_{{z∈Z \atop q(z)>μ}} q(z)r_{z}  (46)
≥ 1 − ∑_{{z∈Z \atop q(z)≤μ}} |π̃_{z} − q(z)| − r̂ ∑_{{z∈Z \atop q(z)>μ}} q(z)
≥ 1 − ϵ − r̂ (using (40) and (43)).  (47)

Observing the Taylor series

\frac{1}{1 − ϵ − r̂} = 1 + ϵ + r̂ + (\frac{1}{1 + θ} )^{3} (ϵ + r̂)^{2} where − ϵ − r̂ ≤ θ ≤ 0,

ρ_{0} ≤ 1 + (1 + \frac{ϵ + r̂}{(1 − ϵ − r̂)^{3}} ) (ϵ + r̂)
≤ 1 + 2.059 … ϵ  (48)

follows.

Vice versa, from (46) one gets

ρ_{0}^{−1} ≤ 1 + ∑_{{z∈Z \atop q(z)≤μ}} |π̃_{z} − q(z)| + r̂ ∑_{{z∈Z \atop μ<q(z)}} q(z)
≤ 1 + ϵ + r̂.  (49)

Observing the Taylor series

\frac{1}{1 + ϵ + r̂} = 1 − ϵ − r̂ + (\frac{1}{1 + θ} )^{3} (ϵ + r̂)^{2} where 0 ≤ θ ≤ ϵ + r̂,

using (49), one obtains

ρ_{0} ≥ 1 − ϵ − r̂ = 1 − 2.016ϵ.

Case where there is a state z′ ∈ Z with q_{z′} > 1 − ϵ: Here, π̃_{z′} = 1.

ρ_{0}^{−1} = ∑_{{z∈Z \atop q(z)≤μ}} π̃_{z} + ∑_{{z∈Z \atop μ<q(z)≤1−ϵ}} q(z)(1 − r_{z}) + 1  (50)
≤ 1 + ϵ + (1 + r̂) ∑_{{z∈Z \atop μ<q(z)≤1−ϵ}} q(z)
≤ 1 + 2ϵ + ϵr̂  (51)

using (39), (23), and ∑_{{z∈Z \atop μ<q(z)≤1−ϵ}} q(z) ≤ 1 − q_{z′} ≤ ϵ. Observing the Taylor series

\frac{1}{1 + 2ϵ + ϵr̂} = 1 − 2ϵ − ϵr̂ + (\frac{1}{1 + θ} )^{3} (2ϵ − ϵr̂)^{2} where 0 ≤ θ ≤ 2ϵ + ϵr̂,

one obtains from (51)

ρ_{0} > 1 − 2.010 … ϵ.  (52)

Vice versa, from (50) one gets

ρ_{0}^{−1} ≥ 1 + ∑_{{z∈Z \atop q(z)≤μ}} π̃_{z} + (1 − r̂) ∑_{{z∈Z \atop μ<q(z)≤1−ϵ}} q(z)
≥ 1.  (53)

In each of the two cases,

1 − 2.016ϵ ≤ ρ_{0} ≤ 1 + 2.059 … ϵ  (54)

holds, and if there is a state z′ ∈ Z with q_{z′} > 1 − ϵ

1 − 2.010 … ϵ < ρ_{0} ≤ 1  (55)

holds.

Now we are prepared to give the approximation errors:

Case q(z) ≤ μ:

π(z) = π̃_{z}ρ_{0} = μρ_{0} ≤ \frac{ϵ}{n}(1 + 2.059 … ϵ)
≤ 1.020 … ϵ/n,

using (25), (39), (54). Clearly,

|π̃_{z} − q(z)| < 1.1ϵ/n

holds.

Case μ < q(z) ≤ 1 − ϵ: In turn, using (24), (41)

\frac{q(z) − π(z)}{q(z)} = \frac{q(z) − π̃_{z}ρ_{0}}{q(z)} = \frac{q(z) − π̃_{z} − π̃_{z}(ρ_{0} − 1)}{q(z)}
= r_{z} − \frac{π̃_{z}}{q(z)}(ρ_{0} − 1) = r_{z} − (1 − r_{z})(ρ_{0} − 1)

holds, and with (43), (54) one gets

\frac{|q(z) − π(z)|}{q(z)} ≤ r̂ + (1 + r̂)2.059 … ϵ
< 3.1ϵ.

Case q(z) > 1 − ϵ: From π(z) = π̃_{z}ρ_{0} = ρ_{0}, and using (55),

|q(z) − π(z)| ≤ 2.1ϵ

follows. This completes the proof of (20).

From (20), (19) follows directly. □

Acknowledgements

The author would like to thank B. Bärk who has been engaged successfully in the DA project for years, W. Everling and M. Heckler for valuable discussions, and the students who worked at the DA idea.

References

[1] B. Bärk, Aggregationsmethoden für Markovketten, 1996.

[2] O. J. Boxma, Workloads and waiting times in single-server systems with multiple customer classes, Queueing Systems 5 (1989) 185–214.

[3] G. Ciardo and K. S. Trivedi, A decomposition approach for stochastic reward net models, Performance Evaluation 18 (1) (1993) 37–59.

[4] J. F. Gierse, Zustandsraum-Aggregation für fast-unabhängige stochastische Petrinetze, Master’s thesis, Universität Bonn, 1995.

[5] T. Heinrichs, B. Bärk, and J. C. Strelen, Wartezeiten für Pollingsysteme mittels numerischer Modelle, in: B. Walke and O. Spaniol, eds., Messung, Modellierung und Bewertung von Rechen- und Kommunikationssystemen (Springer, Berlin, 1993) 173–185.

[6] A. Jennings and W. J. Stewart, Simultaneous iteration for partial eigensolution of real matrices, J. Inst. Maths Applics 15 (1975) 351–361.

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

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

[9] R. Nitzgen, Aggregations-Disaggregationsverfahren für verallgemeinerte stochastische Petrinetze, Master’s thesis, Universität Bonn, 1993.

[10] J. M. Ortega and W. C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables (Academic Press, New York, London, 1970).

[11] B. Philippe, Y. Saad, and W. J. Stewart, Numerical methods in Markov chain modeling, Operations Reseach 40 (6) (1992).

[12] P. J. Schweitzer, A survey of aggregation-disaggregation in large Markov chains, in: W. J. Stewart, ed., Numerical Solution of Markov Chains (Marcel Dekker, New York, Basel, Hong Kong, 1991) 63–88.

[13] J. E. Shore and R. W. Johnson, Axiomatic derivation of the principle of maximum entropy and the principle of minimum cross-entropy, IEEE Trans. Inform. Th. 26 (1) (1980) 26–37.

[14] W. J. Stewart, A comparison of numerical techniques in Markov modeling, Comm. ACM 21 (2) (1978) 144–152.

[15] J. C. Strelen, A fast simultaneous iteration technique for the analysis of Markov chains, Interner Bericht II/89/1, Institut für Informatik, Universität Bonn, 1989.

[16] J. C. Strelen, Iterative Analyse von Markov-Modellen mit alternierender Aggregation und Disaggregation, in: A. Lehmann and F. Lehmann, eds., Messung, Modellierung und Bewertung von Rechensystemen IFB 286 (Springer, Berlin, 1991) 320–336.

[17] J. C. Strelen and B. Bärk, An approach to the numerical analysis of multiple-queue cyclic service systems, in: G. Stiege and J. S. Lie, eds., Messung, Modellierung und Bewertung von Rechensystemen und Netzen IFB 218 (Springer, Berlin, 1989) 75–88.

[18] J. C. Strelen, B. Bärk, J. Becker, and V. Jonas, Analysis of queueing networks with blocking using a new aggregation technique, Submitted for publication.