2. PRELIMINARIES ON THE STOCHASTIC COMPARISON

Here we state the basic definitions and theorems on the stochastic comparison approach. We refer to [13,16] for further details. First, we give the definition to compare two random variables taking values on a totally ordered space [Escr ]. Let [Fscr ]_st denote the class of all increasing real functions on [Escr ] and [Fscr ]_icx denote the class of all increasing and convex real functions on [Escr ]. We denote by ≤_{[Fscr ]} the stochastic order relation, where [Fscr ] can be replaced by st or icx to be associated respectively to the class of functions [Fscr ]_st or [Fscr ]_icx. In the sequel, all of the vectors and matrices are boldfaced, all of the vectors are row vectors, and x^t denotes the transposed vector for x. When comparing two vectors or matrices, ≤ denotes the componentwise comparison.

Definition 1: Let X and Y be two random variables taking values on a totally ordered space [Escr ],

whenever the expectations exist.

Notice that for discrete random variables X and Y with probability vectors p and q, the notations p ≤_{[Fscr ]} q and X ≤_{[Fscr ]} Y are used interchangeably. Stochastic orders ≤_st and ≤_icx can be also defined by matrices (see [9,10]). Throughout this work, we assume that [Escr ] = {1,…, n}, but the following statements can be extended to the infinite case. We give K_st and K_icx matrices related respectively to the ≤_st and ≤_icx orders. In the sequel, K_{[Fscr ]} denotes the matrix related to the ≤_{[Fscr ]} order, [Fscr ] ∈ {st,icx}:

Proposition 1: Let X and Y be two random variables taking values on [Escr ] = {1,2,…, n} and let p = (p_i)_i=1ⁿ and q = (q_i)_i=1ⁿ, be probability vectors such that

which can be also given as follows:

It is shown in Theorem 5.2.11 of [13, p.186] that monotonicity and comparability of the probability transition matrices of time-homogeneous Markov chains yield sufficient conditions to compare the underlying chains stochastically. We first define the monotonicity and comparability of stochastic matrices and then state this theorem.

Definition 2: Let P be a stochastic matrix. P is said to be stochastically ≤_{[Fscr ]}-monotone if for any probability vectors p and q,

Definition 3: Let P and Q be two stochastic matrices. Q is said to be an upper-bounding matrix of P in the sense of the ≤_{[Fscr ]} order (P ≤_{[Fscr ]} Q) if

Let us remark that this is equivalent to saying that P ≤_{[Fscr ]} Q if

where P_i,* denotes the ith row of matrix P.

Theorem 1: Let P (resp. Q) be the probability transition matrix of the time-homogeneous Markov chain {X_n,n ≥ 0} (resp. {Y_n,n ≥ 0}). If

• X₀ ≤_{[Fscr ]} Y₀,
• at least one of the probability transition matrices is monotone (i.e., either P or Q is ≤_{[Fscr ]}-monotone),
• the transition matrices are comparable (i.e., P ≤_{[Fscr ]} Q),

then

The following lemma [13] lets us compare the steady-state distributions of Markov chains, which can be considered as the limiting case (Π_P = lim_n→∞ X(n)).

Lemma 1: Let Q be a monotone, upper-bounding matrix for a finite stochastic matrix P in the sense of the ≤_{[Fscr ]} order. If the steady-state distributions (π_P, π_Q) exist, then

3. GENERALIZED CLASS [Cscr ] ([Cscr ]^{[Gscr ]}) MARKOV CHAINS

In this section we give the generalization of class [Cscr ] Markov chains that will be called class [Cscr ]^{[Gscr ]} Markov chains. First we present definitions and main properties of this class; then we give algorithms to construct bounding chains.

3.1. Class [Cscr ]^{[Gscr ]} Markov Chains and Closed-Form Solutions

Definition 4 (Class [Cscr ]^{[Gscr ]}): A stochastic matrix P belongs to [Cscr ]^{[Gscr ]} if there are three vectors v, c, and r in

such that v is a probability vector (v_i ≥ 0, ∀i,

,

, and

where e = (e_i)_i=1ⁿ and e_i = 1, ∀i.

Remark 1: Class [Cscr ] Markov chains [3] can be obtained by taking r = (r_i)_i=1ⁿ, r_i = i − 1, ∀i. In this case, v equals P_1,*.

In the following lemma we give the properties of power matrices for class [Cscr ]^{[Gscr ]} matrices.

Lemma 2: Let P be a stochastic matrix from class [Cscr ]^{[Gscr ]}: P = e^tv + r^tc. Let α and β be the following constants:

The power matrices P^k,k ≥ 2, belong to class [Cscr ]^{[Gscr ]} and they can easily be computed as follows:

Proof: The proof is done by induction on k.

Basic step (k = 2): P² = (e^tv + r^tc)² = e^tve^tv + e^tvr^tc + r^tce^tv + r^tcr^tc. Using the fact that ve^t = 1 and ce^t = 0 (Definition 4), cr^t = α and vr^t = β (Eq. (4)), we have: P² = e^tv + [αr^t + βe^t] c.

Induction step: Suppose that the lemma holds for k and let us show that it holds also for k + 1. Let

. Hence, P^k = e^tv + f^tc, so

On the other hand,

and

. Therefore,

. █

We can now state the theorem on the closed-form solution for transient distributions of this class of Markov chains.

Theorem 2: Let {X(k), k ≥ 0} be a time-homogeneous discrete-time Markov chain with probability transition matrix P and let π^k = (π_i^k)_i=1ⁿ be the distribution of X(k).

If P belongs to class [Cscr ]^{[Gscr ]}, then for all k ≥ 0,

where a_k is the constant defined as

The constants α and β have been already defined in Eq. (4) and g is defined as

Proof: Using Lemma 2, we have:

The following corollary gives the conditions to have the limiting distribution of [Cscr ]^{[Gscr ]} matrices.

Corollary 1 (Steady-state distribution for class [Cscr ]^{[Gscr ]} matrices): Let {X(k), k ≥ 0} be a time-homogeneous discrete-time Markov chain with probability transition matrix P from class [Cscr ]^{[Gscr ]} and let α be defined as α = cr^t. {X(k), k ≥ 0} has a stationary distribution if and only if |α| < 1. In that case, the unique stationary distribution (steady-state distribution) π = lim_k→∞ π^k is given by

where R is a constant given by R = vr^t/(1 − cr^t).

Proof: From Theorem 2 we have π^k = v + a_k c, ∀k, where π^k is the distribution of X(k) and

. For |α| < 1, lim_k→∞ a_k = β/(1 − α), so we have π = v + (β/(1 − α))c = v + Rc, which is independent of π⁰. For |α| ≥ 1, a_k diverges; thus, the stationary distribution does not exist. █

Let us now give the conditions to have the limiting distribution in terms of entries of P belonging to class [Cscr ]^{[Gscr ]}.

Remark 2: Since

, we have α ≥ −1. The condition |α| < 1 is satisfied if and only if

.

It follows from Theorem 2 and Corollary 1 that the computation of transient distributions and the steady-state distribution are linear with the size of the model.

Remark 3: The complexity of computing π^k is θ(max{n,k}) and it equals θ(n) for π.

3.2. Monotonicity of Class [Cscr ]^{[Gscr ]} Matrices

The following proposition gives sufficient conditions for the monotonicity of class [Cscr ]^{[Gscr ]} matrices in terms of vectors c and r. For a vector

, we will write x ∈ [Fscr ] if and only if vector x, seen as a real function on [Escr ] = {1,…, n}, is in class [Fscr ].

Proposition 2: A matrix P = e^tv + r^tc ∈ [Cscr ]^{[Gscr ]} is ≤_{[Fscr ]}-monotone if

Proof: First, we define vectors c′ = cK_{[Fscr ]} and v′ = vK_{[Fscr ]}. Hence, PK_{[Fscr ]} = e^tv′ + r^tc′. To prove that P is ≤_{[Fscr ]}-monotone, we consider two probability vectors p and q such that p ≤_{[Fscr ]} q. According to Definition 2 we must show that pP ≤_{[Fscr ]} qP, which is equivalent to pPK_{[Fscr ]} ≤ qPK_{[Fscr ]} (see Proposition 1). We have

On the other hand, if X and Y denote random variables corresponding respectively to probability vectors p and q, then using Definition 1 and the fact that r ∈ [Fscr ], we have E [r(X)] ≤ E [r(Y)]. Remark that

and, similarly, E [r(Y)] = qr^t, so pr^t ≤ qr^t. Using in addition the fact that c′ = cK_{[Fscr ]} ≥ 0, we conclude that pr^tc′ ≤ qr^tc′, and the proof is achieved. █

Remark 4: For the ≤_st order, r ∈ [Fscr ]_st if and only if vector r is increasing; that is,

For the ≤_icx order, r ∈ [Fscr ]_icx if and only if vector r is increasing and convex; that is,

In the general case, when the underlying matrix P does not belong to class [Cscr ]^{[Gscr ]}, the monotonicity characterization is given as follows:

• P is ≤_st-monotone if and only if for each column j, the function
  
  is increasing with respect to i [9].
• P is ≤_icx-monotone if and only if for each column j, the function
  
  is increasing and convex with respect to i (see [10] for the infinite state space case and [2] for the finite state space case).

3.3. Bounding Algorithms

Obviously, having closed-form solutions makes [Cscr ]^{[Gscr ]} matrices attractive for computational issues. However, for a given matrix, the constraints to belong to [Cscr ]^{[Gscr ]} matrices might not be satisfied. We propose to construct class [Cscr ]^{[Gscr ]} bounding matrices by applying the stochastic comparison approach. In fact, we apply Theorem 1: For a given P, we construct Q ∈ [Cscr ]^{[Gscr ]} such that P ≤_{[Fscr ]} Q and Q is ≤_{[Fscr ]}-monotone. Therefore, we can compute the transient distributions and the steady-state distribution of Q through its closed-form solutions to provide bounding distributions on P.

3.3.1. Strong stochastic bounds.

In this subsection we give the algorithm to construct bounding matrix in the sense of the ≤_st order. As stated earlier, a [Cscr ]^{[Gscr ]} stochastic matrix is defined through three vectors: v, r, and c. For a given matrix P and a nonnegative, increasing vector r, Algorithm 1 returns vectors c and v such that Q = e^tv + r^tc ∈ [Cscr ]^{[Gscr ]}; P ≤_st Q and Q is ≤_st-monotone. The nonnegativity constraint is used only to simplify the algorithm and its proof. Some possible choices for r will be discussed in Section 3.4.

Algorithm 1: ≤_st-Monotone [Cscr ]^{[Gscr ]} Upper Bound Q of P

Input: P: stochastic matrix; r: nonnegative, r ∈ [Fscr ]_st (increasing)

Output: Q: ≤_st-monotone, [Cscr ]^{[Gscr ]} upper bound of P, Q = e^tv + r^tc

Notation: K := min{i | r_i > 0},M := min{i | r_i = r_n}

We suppose that M ≠ 1. (In the other case, the algorithm is trivial: Taking the maximal value of each column of PK_st, we obtain a vector x. We can then take v = xK_st⁻¹; Q = e^tv.)

For M ≠ 1, K is well defined and we have K ≤ M.

//Column n:

1.

;

2.

;

for j = n − 1 downto 2 do

3.

;

4.

;

end

//Column 1:

5.

;

6.

;

We give the following theorem to prove Algorithm 1.

Theorem 3: Let P be a stochastic matrix and Q be the matrix computed by Algorithm 1. The matrix Q is a stochastic matrix, Q ∈ [Cscr ]^{[Gscr ]}, Q is ≤_st-monotone and P ≤_st Q.

The proof of Theorem 3 is given in the Appendix.

3.3.2. Increasing convex bounds.

In Algorithm 2 we give the construction of ≤_icx bounding matrices. Let us first define the operators that will be used in this algorithm.

and

are defined as follows:

Let us remark that φ_icx(x,j) = (xK_icx)_j and

.

Algorithm 2: ≤_icx-Monotone [Cscr ]^{[Gscr ]} Upper Bound Q of P

Input: P: stochastic matrix; r: nonnegative vector, r ∈ [Fscr ]_icx (increasing, convex)

Output: Q: ≤_icx-monotone [Cscr ]^{[Gscr ]} upper bound of P, Q = e^tv + r^tc

Notation: K := min{i | r_i > 0},M := min{i | r_i = r_n}

We suppose that M ≠ 1. (In the other case, the algorithm is trivial: Taking the maximal value of each column of PK_icx, we obtain a vector x. We can then take v = xK_icx⁻¹; Q = e^tv.)

For M ≠ 1, K is well defined and we have K ≤ M.

//Column n:

1.

;

2.

;

for: j = n − 1 downto 2 do

3.

;

4.

;

end

//Column 1:

5.

;

6.

;

The following theorem is given to prove Algorithm 2.

Theorem 4: Let P be a stochastic matrix and Q be the matrix computed by Algorithm 2. The matrix Q is a stochastic matrix from class [Cscr ]^{[Gscr ]} and it is an ≤_icx-monotone upper bound for P.

The proof of Theorem 4 is given in the Appendix.

Remark 5: The computation complexity of Algorithm 1 and Algorithm 2 is θ(n²).

3.4. Choosing Vector r

For both algorithms, vector r of row increments, satisfying monotonicity conditions (see Proposition 2), is fixed. Moreover, vector r is also supposed to be nonnegative in order to simplify the algorithm and its proof. Thus, any vector satisfying the above conditions can be chosen as the vector of row increments. We present some of possible choices for vector r in the following.

If we take vector r = (r_i)_i=1ⁿ, r_i = i − 1, we obtain the bounding matrix that belongs to class [Cscr ] (see Remark 1). It is important to notice that this choice is independent of the underlying matrix entries.

We propose some other possibilities for vector r in order to compute more accurate bounds. The main idea here is to improve bounding entries in the last column. In fact, since in both algorithms the construction is done recursively from the last column to the first column, the perturbations in the last column entries of the bounding matrix are propagated to other bounding matrix entries. So by diminishing the perturbations in the last column, it might be possible to improve the accuracy of the bounds. The other reason comes essentially from performance evaluation studies. Quite often, we are interested in bounding only the tail distribution. For instance, in reliability applications we are interested in bounding the probability that the studied system is not operational; in queuing applications we bound the overflow probability of the buffer. In general, these measures concern a few states and they can be ordered as last states. Thus, accurate bounds for last states would be sufficient for these cases.

3.4.1. Strong stochastic case.

For the strong stochastic order, there is a vector r that minimizes the differences between the exact and the bounding entries for the last column.

Proposition 3: Let the input vector r of Algorithm 1 be defined as

Then the bounding matrix Q computed by Algorithm 1 satisfies the following:

Proof: P ≤_st U and P ≤_st Q (see Theorem 3) yield that p_i,n ≤ u_i,n and p_i,n ≤ q_i,n, ∀i, so Eq. (9) is equivalent to q_i,n ≤ u_i,n, ∀i. Since U is ≤_st-monotone, u_i,n are increasing in i. Moreover, u_i,n ≥ p_i,n, ∀i, so u_i,n ≥ max_k≤i p_k,n = r_i, ∀i.

Finally, we will show that Algorithm 1 with the input vector r defined by Eq. (9) yields v_n = 0 and c_n = 1. Therefore, q_i,n = v_n + r_i c_n = r_i ≤ u_i,n, ∀i, and the proof is concluded. Indeed, we have r_n − r_i > 0, 1 ≤ i < M, r_i = max_k≤i p_k,i ≥ p_i,n, ∀i, and r_n ≤ 1; thus, it can be seen from the first line of Algorithm 1 that v_n = 0. On the other hand, p_i,n ≤ r_i, ∀i, and r_K = p_K,n so it follows from the second line of Algorithm 1 that c_n = 1. █

In the sequel, vector r obtained by Eq. (8) will be called the St-LC vector. In the following example, we present class [Cscr ]^{[Gscr ]} matrices obtained with vector r corresponding to class [Cscr ] and with the St-LC vector computed from Eq. (8).

Example 1: Let P be the original matrix with the steady-state distribution

. By taking vector r that corresponds to class [Cscr ], r = (0, 1, 2), Algorithm 1 returns matrix A with vector v = (0.4,0.6,0) and vector c = (−0.2, −0.2, 0.4). The steady-state distribution of A is π_A = (0.1, 0.3, 0.6).

The St-LC vector computed from Eq. (8) is r = (0, 0.4, 0.6). In this case, Algorithm 1 returns bounding matrix A′ with vector v = (0, 0.4, 0.6) and vector c = (−0.5, −0.5, 1). The steady-state distribution of this bounding matrix is π_A′ = (0.2, 0.4, 0.4).

We note that A′ belongs to class [Cscr ]^{[Gscr ]} but not to class [Cscr ] and π_P ≤_st π_A′ ≤_st π_A. However, we do not always have this relation (see Fig. 4).

Steady-state analysis (≤st order).

3.4.2. Increasing convex case.

In the case of the increasing convex order, for a given matrix P it is not always possible to find a monotone upper-bounding matrix Q satisfying the following conditions:

Let us show this for matrix P of Example 1. It can be shown that the following matrices B and B′ are ≤_icx-monotone upper-bounding matrices for P. To satisfy Eq. (10) we must have p_i,n ≤ q_i,n ≤ min(b_i,n,b_i,n′), ∀i; thus, q_1,3 = 0, q_2,3 = 0.4, and q_3,3 = 0.6. However, this vector r = (0, 0.4, 0.6) is not compatible with the monotonicity constraints (see Remark 4).

Let us recall that an ≤_icx-monotone upper-bounding matrix Q of P is optimal in the sense of ≤_icx order if

In fact, Eq. (10) represents the optimality conditions for the last column, so it is a necessary condition for the global optimality [Eq. (11)]. Hence, the previous example shows that it is not always possible to find the optimal bound in the ≤_icx sense. Moreover, this is the case for any stochastic matrix, not only for class [Cscr ]^{[Gscr ]} matrices, since the properties of this class have not been considered.

As we cannot always satisfy Eq. (10), we propose to minimize the sum of distances between q_i,n and p_i,n:

Proposition 4: Let vector r be a solution of the following problem: Find vector x with a minimal sum of entries satisfying the monotonicity (x ∈ [Fscr ]_icx) and the comparison constraints,

Then the upper-bounding matrix Q computed by Algorithm 2 with this input vector r minimizes Eq. (12); that is,

Proof: Let us remark first that P ≤_icx U and P ≤_icx Q (see Theorem 4) give p_i,n ≤ u_i,n and p_i,n ≤ q_i,n, ∀i; so Eq. (14) is equivalent to

. Constraints (13) are fulfilled by all ≤_icx-monotone upper bounds of P, so we have

, ∀U ≤_icx-monotone upper bound of P.

Similarly to the proof of Proposition 3, it can be shown from Algorithm 2 that v_n = 0 and c_n ≤ 1. Moreover, if c_n < 1, then

, which is in contradiction with the fact that r is a solution of (13). Thus, we have c_n = 1. █

The Simplex algorithm can be used to find a solution of problem (13). Let us remark that solution of this problem is not necessarily unique. The bounds computed from Algorithm 2 with input vector r is obtained through the Simplex algorithm will be called Icx-Simplex bounds.

Although problem (13) can be efficiently solved by the Simplex algorithm, we cannot guarantee its complexity in the worst case. Therefore, we present an alternative heuristic choice for r. Unrolling the comparison [Eq. (2)] and the monotonicity (Remark 4) constraints for the last column increasingly in rows,

we obtain a nonnegative vector r such that r ∈ [Fscr ]_icx. In the sequel, the bounds obtained from Algorithm 2 with this input vector will be called Icx-LC bounds.

Let us remark that a vector r obtained by relations (15) satisfy constraints (13) except the last one, r_n ≤ 1. However, this will be compensated by the fact that we will have c_n < 1 during the construction of the bounding matrix through Algorithm 2. Thus, the last column of the bounding matrix Q computed by Algorithm 2 satisfies all of the constraints (13). Although the Icx-LC bound does not necessarily minimize Eq. (12), since we take into consideration the information contained in the last column of the original matrix, this choice could yield better bounds than class [Cscr ] bounds for rewards that are mostly based on the last states. This is true in the case of the example studied in Section 4. Moreover, for this example, Icx-LC bounds are remarkably close to Icx-Simplex bounds.

Remark 6: The complexity of computing different vectors r is not important compared to that of Algorithms 1 and 2 (Remark 5). In fact, Eqs. (8) and (15) can be solved with linear complexity, and Eq. (13) can be solved efficiently with the Simplex algorithm.

4. NUMERICAL EXAMPLE

Let us first emphasize that the proposed bounding method can be applied to any finite discrete-time Markov chain, but the tightness of bounds depends largely on the underlying transition matrices. Our main goal in this section is to give some insights on the relationship between the accuracy of [Cscr ]^{[Gscr ]} bounds and the underlying matrix structures. We study an example that has been also studied in previous articles because of its simplicity in generating different structures of probability transition matrices.

We consider a finite-capacity buffer with multiple servers and batch arrivals. The time is supposed to be divided into slots and the analysis is provided in discrete time. The number of servers is denoted by S and the buffer capacity by B; thus, C = B + S is the maximum number of customers in the system (see Fig. 1). The size of batch arrivals is assumed to be distributed according to a truncated modified geometrical distribution with parameter p_a. The probability that i customers arrive during a slot is a_i = (1 − p_a)p_aⁱ, 1 ≤ i ≤ A, where A is the maximal batch size. The service time is distributed geometrically, and p_s denotes the probability that the service is finished in one slot. Obviously, the values of S and A have a large impact on the sparseness of the matrix. The number of nonzero entries increases with the values of S and A and the probability transition matrix is full when A = S = B. The probabilities p_a and p_s for arrivals and services, respectively, determine where the probability mass is located in the matrix. In the following figures the curves of bounds are noted according to vector r used as input in bounding algorithms (see Section 3.4).

Considered system.

In Figure 2, the values of B, p_a, and p_s are fixed and we assume that S = A. We vary this last value and give the average number of customers in the steady state. We note that when S = A increases, the accuracy of bounds is improved for all class [Cscr ]^{[Gscr ]} bounds for both ≤_st and ≤_icx orders. In fact, whatever the structure of the original matrix, the bounding matrix is filled due to the structure of class [Cscr ]^{[Gscr ]} matrices. Thus, if the number of zero entries in the original matrix is large, the number of modified entries would be also large and this might cause a deterioration of the accuracy of bounds. However, we note that it is possible to have tight bounds even for relatively sparse matrices, especially with the ≤_icx order. In fact, the presented results confirm the intuition which follows from the structure of [Cscr ]^{[Gscr ]} matrices that the proposed bounds are more adapted for full matrices. On the other hand, it is not possible to state that the bounds for relatively sparse matrices are always loose, since the quality of bounds depends also on the magnitude of the perturbations.

Impact of sparseness of the transition matrix on the quality of class [Cscr ][Gscr ] bounds (steady state).

In the case of the ≤_st order, we can see in Figure 2 that the St-LC and the St-Class C curves cross each other: The class [Cscr ] bound is tighter for large values of A = S (full matrices) and the St-LC bound is tighter for sparser matrices. Thus, each of these ≤_st bounds might be of interest. In the case of the ≤_icx order, class [Cscr ] bounds are tighter than the other considered class [Cscr ]^{[Gscr ]} bounds and they are very close to the exact values for full matrices. Moreover, there is not a sensible difference between Icx-LC and Icx-Simplex bounds. This phenomena has been observed not only in Figure 2 but also for different sets of parameters that we considered. Hence, the heuristic proposed in Section 3.4.2 seems to be useful to avoid the complexity of the Simplex algorithm.

We also note in Figure 2 that ≤_icx bounds are generally more accurate than the ≤_st bounds. This is due to specific monotonicity properties of class [Cscr ]^{[Gscr ]} matrices. In fact, the ≤_st comparison of random variables and stochastic matrices always implies the ≤_icx comparison [13], as the class of increasing convex functions is contained in the class of increasing functions. For the monotonicity, this is not generally the case. However, for class [Cscr ] matrices, the monotonicity in the ≤_st sense implies the monotonicity in the ≤_icx sense [4]. We can also show that, for [Cscr ]^{[Gscr ]} matrices, if the chosen vector r is increasing and convex (r ∈ [Fscr ]_icx), then cK_st ≥ 0 implies the ≤_icx monotonicity (see Proposition 2). On the other hand, ≤_st bounds might be useful if one wants to compute bounds on increasing but not convex rewards. For instance, if we want to compare the tail distributions, we cannot use ≤_icx bounds.

Let us remember that the proposed approach can be also applied to provide transient bounds. In Figure 3 we present transient behaviors of the average number of customers for the parameter values of Figure 2 with A = S = 10. We can see here that in the transient case we can make the same observations as in the steady-state case (Figure 2). More generally, we have noticed also through other experiments that the quality of transient bounds has the same tendency as in the case of steady-state bounds.

Transient [Cscr ][Gscr ] bounds for the case A = S = 10.

In Figures 4 and 5, we consider bounds on steady-state distributions. Let N be the number of customers in the steady-state case. In Figure 4, we present ≤_st bounds for the tail distribution (Prob(N ≥ x)). For ≤_icx bounds, it is not possible to compare tail distributions. However, we have the following characterization of the ≤_icx order (see [13]):

Steady-state analysis (≤icx order).

Informally speaking, this characterization says that, for a given level x, we compare the mean overflow from that level. In Figure 5, we compare E [(N − x + 1)⁺] for different ≤_icx bounds. We note that for this set of parameters, the bounds obtained using vectors r derived from the last column are more accurate for both ≤_icx and ≤_st bounds. For instance, in Figures 4 and 5, for x ≥ 10 we obtain tighter bounds by taking into account the last column of the original matrix.

It is not surprising to obtain more accurate bounds for last states (i.e., bounds on the rewards of last states), by taking into account the information contained in the last column of the original matrix. In fact the information used to derive the input vector r is valid for the last column and also, to a lesser extent, for the columns that are relatively close to it. However, it is less intuitive that we obtain tighter results from class [Cscr ] ≤_icx bounds when studying average behavior of the underlying model.

It can be seen from Figures 6 and 7 that different class [Cscr ]^{[Gscr ]} bounds can cross each other several times. Due to the low computational complexity of class [Cscr ]^{[Gscr ]} bounds, it is worth computing different class [Cscr ]^{[Gscr ]} bounds and taking the best one for each value. In Table 1 we give the best ≤_st and ≤_icx bounds that can be derived from the bounds of Figures 6 and 7. In Table 2 we present, for the same set of parameters, the best transient bounds for Prob(N ≥ 15) for the ≤_st order and E [(N − 14)⁺] for the ≤_icx order. We observe that the best transient bounds are obtained from different vectors r, as in the steady-state case.

Bounds can cross more than once (≤st case).

Bounds can cross more than once (≤icx case).

Best ≤st and ≤icx Bounds Derived from Bounds of Figures 6 and 7

Best ≤st and ≤icx Transient Bounds for the Same Set of Parameters as in Figures 6 and 7

5. CONCLUSION

In this article we present the generalization of class [Cscr ] matrices to class [Cscr ]^{[Gscr ]} matrices. This class of Markov chains has closed-form solutions to compute transient distributions and the steady-state distribution, which makes it attractive for computational issues of large Markov chains. We propose to apply the algorithmic stochastic comparison approach to construct upper-bounding matrices belonging to this class. For a given finite Markov chain we construct a bounding class [Cscr ]^{[Gscr ]} chain and compute its distributions to provide bounding distributions. By analyzing bounding class [Cscr ]^{[Gscr ]} chains, we significantly gain in terms of both computational and storage complexity, compared to the exact analysis.

We give algorithms to construct class [Cscr ]^{[Gscr ]}, monotone, upper-bounding matrices in the sense of the ≤_st and the ≤_icx orders. Class [Cscr ]^{[Gscr ]} matrices are completely defined by means of three vectors: v, c, and r. Our method consists of two steps. First, we use the information contained in the last column of the original transition matrix to compute vector r of row increments. Once vector r is fixed, we apply Algorithm 1 (≤_st bounds) or Algorithm 2 (≤_icx bounds) to determine the remaining v and c vectors.

We present different choices of vector r to construct class [Cscr ]^{[Gscr ]} bounding matrices in the sense of both ≤_st and ≤_icx orders. By changing input vector r, it is possible to easily construct different bounding matrices in class [Cscr ]^{[Gscr ]}. These matrices can be analyzed through the closed-form solutions and the better bound can be taken for a given performance measure of interest.

We have performed numerical experiments to study the accuracy of bounds as a function of matrix structures and different choices of vector r. We observe that class [Cscr ]^{[Gscr ]} bounds become tighter when the sparseness of the underlying matrix decreases and that ≤_icx bounds are tighter in general than ≤_st bounds. We show that ≤_st bounds can be improved by class [Cscr ]^{[Gscr ]} matrices compared to class [Cscr ] matrices. This is especially useful when we are interested in functionals of distributions that are increasing but not convex. For ≤_icx bounds, we conclude that class [Cscr ] bounds are generally accurate. Other class [Cscr ]^{[Gscr ]} bounds could be interesting when computing bounds for rewards defined on last states.

Computing closed-form stochastic bounds seems to be a promising approach to overcome the state space explosion problem of Markov chains. We are considering extending our work to look for other matrix structures having closed-form solutions and to other stochastic orders.